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1.  Introduction: 


Optoacoustic  tomography  (OAT),  also  known  as  photoacoustic  tomography,  is  a  rapidly 
emerging  hybrid  imaging  modality  that  combines  high  ultrasonic  resolution  with  strong  optical 
or  microwave  contrast,  and  represents  a  highly  promising  biomedical  imaging  modality  with  a 
great  potential  for  early  breast  cancer  detection  and  subsequent  management.  On  the  other  hand, 
ultrasound  computed  tomography  (UST)  can  provide  high-resolution  anatomical  images  of 
breast  lesions  based  on  three  complementary  acoustic  properties  (speed-of-sound  (SOS), 
attenuation,  reflectivity).  The  broad  objective  of  the  proposed  research  is  to  develop  an 
optimized  system  design  and  associated  image  reconstruction  algorithms  for  a  hybrid  three- 
dimensional  (3D)  breast  imaging  system  that  combines  OAT  and  UST.  This  system  will  provide 
co-registered  functional  and  anatomical  information  without  the  need  for  breast  compression, 
ionizing  radiation,  or  contrast  agents,  which  will  result  in  a  highly  effective  and  affordable  3D 
imaging  method. 

2.  Keywords: 

Ultrasound  Computed  tomography,  Optoacoutics  Tomography,  Breast  Imaging,  Iterative  Image 
reconstruction  algorithms,  Bent-Ray  USCT,  GPU-accelerated  algorithms 

3.  Overall  Project  Summary 

Task  1:  Construct  a  computational  model  of  the  OAT/UST  imager  and  identify  optimal 
system  geometries: 

la.  Computational  modeling  of  the  proposed  imager:  We  have  conducted  a  computational 
study  of  imager  designs  for  3D  laser  ultrasound  (LU)  UST.  The  imager  consists  of  a  circular 
transducer  arc  attached  opposite  to  a  half  cylinder.  There  were  128  transducers  attached  to  the 
circular  arc;  which  has  140  mm  diameter  and  118  mm  aperture  size  and  is  symmetric  around  the 
central  plane. 


lb.  Optimization  studies  (Months  8-36):  The  goal  of  the  optimization  studies  was  to  identify 
optimal  imaging  system  configurations  for  OAT  and  UST  that  satisfy  the  design  goals.  The 
number  and  location  of  the  ultrasound  emitters,  which  are  constrained  to  reside  on  this 
cylindrical  surface,  were  optimized.  In  addition  to  the  system  parameters,  an  iterative  image 
reconstruction  was  optimized.  We  demonstrated  that  high  quality  3D  SOS  maps  can  be 
reconstructed  when  only  32  emitters  and  128  receiving  transducers  are  employed  to  record  time- 
of-flight  data  at  360  tomographic  view  angles. 


lc.  Approaches  to  evaluation  of  image  quality  for  system  optimization:  Standard  quantitative 
measures  of  image  quality  were  employed.  These  included  bias-variance  curves,  contrast-to- 
noise  ratios. 

It  is  advocated  in  the  modem  imaging  science  literature  to  utilize  objective,  or  task-based, 
measures  of  system  performance  to  guide  the  optimization  of  hardware  design  and  image 
reconstmction  algorithms.  We  investigated  this  approach  to  assess  the  performance  of  OAT 
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breast  imaging  systems.  In  particular,  we  applied  principles  from  signal  detection  theory  to 
compute  the  detectability  of  a  simulated  tumor  at  different  depths  within  a  breast,  for  two 
different  system  designs.  The  signal-to-noise  ratio  of  the  test  statistic  computed  by  a  numerical 
observer  was  employed  as  the  task-specific  summary  measure  of  system  performance.  A 
numerical  breast  model  was  employed  that  contains  both  slowly  varying  background  and  vessel 
structures  as  the  background  model,  and  superimpose  a  deterministic  signal  to  emulate  a  tumor. 
This  study  quantified  how  signal  detection  performance  of  a  numerical  observer  will  vary  as  a 
function  of  signal  depth  and  imaging  system  characteristics.  The  described  methodology  is  being 
employed  to  systematically  optimize  our  OAT  imaging  system  design  for  tumor  detection  tasks. 


Task  2:  Development  of  reconstruction  methods  for  sparse-array  3D  UST 
Reconstruction  of  SOS  distribution: 

Iterative  image  reconstruction  algorithms  have  been  developed  to  reconstruct  the  SOS 
distribution.  These  algorithms  have  also  been  validated  for  experimental  phantom  studies.  To 
perform  USCT,  measured  time-of-flight  (TOF)  was  extracted  from  the  measured  signals  of  the 
transducer  elements. 

(i)  Time-of-flight  (TOF)  extraction  of  the  transmission  ultrasound  signals:  I  utilized 
geometrical  acoustic-based  ray  theory  to  establish  a  non-linear  model  that  relates  the  measured 
TOF  values  to  speed  of  sound  (SOS)  distribution.  To  solve  this  nonlinear  optimization  problem, 
we  needed  to  extract  time-of-flight  (TOF)  from  the  measured  signal.  This  is  a  very  important 
pre-processing  step  for  good  image  quality.  In  search  of  the  best  TOF  extraction  technique,  six 
TOF  extraction  algorithms  have  been  implemented  and  compared:  (i)  envelope-detection 
method,  (ii)  picking  the  max  value  of  filtered  signal,  (iii)  AlC-method,  (iv)  weighted-AIC  picker, 
(v)  cross-correlation  method,  and  (vi)  thresholding  method  from  windowed  and  filtered  signal 
[1].  These  methods  were  investigated  for  both  the  computer  simulation  by  adding  different  levels 
of  noise  and  the  experimental  data  (provided  by  Tomowave  Laboratories  Inc.,  Houston  TX). 

(ii)  Bent-Ray  method:  Algorithms  are  developed  for  reconstructing  the  SOS  distribution  of 
breast  from  knowledge  of  time-of-flight  (TOF)  measurements  of  the  transmission  ultrasound 
signals.  Utilizing  the  geometrical  acoustic-based  ray  theory,  a  non-linear  model  has  been 
established  that  relates  the  measured  TOF  values  to  the  SOS  distribution. 

For  a  given  SOS  distribution,  numerically  solution  of  the  Eikonal  equation  yielded  the  ray  paths. 
An  iterative  reconstruction  method  was  developed  for  inverting  the  resulting  system  of  equations 
that  alternatively  updates  the  estimates  of  the  SOS  and  ray  paths,  minimizes  a  regularized  cost 
function  to  obtain  the  final  estimate  of  the  SOS  [2], 

(iii)  Adjoint-State  Method:  I  also  investigated  a  partial  differential  equation-based  Eulerian 
approach  to  travel-time  tomography  as  an  alternative  approach  [3].  The  work  on  comparison  of 
the  Adjoint-State  Method  are  Ray-tracing  algorithm  was  presented  in  SPIE  Photonics  West, 
2014  and  SPIE  Medical  Imaging  2014.  For  detail  implementation  of  this  algorithm  please  see 
attached  proceedings  paper.  The  waveform  inversion  method  for  SOS  reconstruction  has  also 
been  explored  in  the  group.  In  this  regard,  adjoint-state  method  provides  a  suitable  initial  SOS 
distribution  to  aid  waveform  inversion  method. 

I  performed  several  numerical  studies  to  compare  bent-ray  and  adjoint-state  method.  The 
method  was  also  validated  for  the  experimental  phantom  studies.  In  this  study,  the  pressure 
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signals  were  recorded  using  a  transducer  array  consisting  of  64  elements  and  a  single  Laser 
Induced  Source  placed  opposite  to  the  middle  transducer.  The  data  was  recoded  for  150  views. 
The  experimental  phantom  consisted  of  three  tubes  at  varying  salt  concentration  to  exhibit 
different  acoustic  and  optical  properties. 


2b.  Reconstruction  of  attenuation  distribution 

Not  accomplished 


2c.  Reconstruction  of  reflectivity:  An  algorithm  has  been  developed  to  produce  reflectivity 
maps.  These  algorithms  are  based  on  the  Synthetic  Transmit  Aperture  (STA)  approach  [4],  This 
method  utilizes  multiple  elements  or  a  single  source  to  produce  the  spherical  waves  and  the 
whole  image  is  being  reconstructed  for  each  emitted  signal.  The  final  reflectivity  map  is  obtained 
by  accumulating  these  individual  images.  It  has  been  shown  that  the  STA  method  improves 
SNR.  This  method  is  especially  useful  in  our  current  study  because  we  are  using  laser  induced 
ultrasound  emitter  (LUS),  which  produce  spherical  waves.  The  reflectivity  map  for  breast 
tomographic  was  not  accomplished. 


Task  3:  Ultrasound-assisted  OAT  image  reconstruction: 

3a.  (i)  Development  of  imaging  models  and  reconstruction  algorithms:  An  interpolation- 
based  discrete-discrete  imaging  model  has  been  implemented  to  perform  3D  OAT  for  breast 
imaging  [5],  In  the  new  implementation,  an  unmatched  back  projection  (or  pixel-driven)  scheme 
has  been  used  and  validated  in  computer  simulations  studies.  This  algorithm  is  five  times  more 
efficient  than  the  ray-driven  back-projection  and  allows  to  perform  iterative  image  reconstruction 
for  large  fields-of-views,  making  it  very  suitable  for  breast  imaging.  To  efficiently  mitigate  data 
incompleteness,  noise,  and  model  error,  I  investigated  the  least-squares  objective  regularized  by 
a  TV-norm  penalty.  I  implemented  the  fast  iterative  shrinkage/thresholding  algorithm  (FISTA)  to 
minimize  cost  function  with  TV  regularization  [6], 

(ii)  In-vivo  OAT-Imaging  to  validate  unmatched  imaging  model:  Numerical  phantom  study 


Figure  1:  Schematic  of  full  body 
mouse  imaging  module.  TomoWave 
Inc.  performed  experiment. 


Figure  2:  Maximum  intensity  projection  renderings  of  3D  images  of 
a  mouse  body  using  (a)  the  original  interpolation-based  algorithm, 
and  (b)  the  accelerated  interpolation-based  algorithm. 


of  the  comparison  between  matched  and  unmatched  imaging  model  was  reported  previously.  It  is 
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important  to  validate  the  model  for  noisy  and  experimental  data.  In  order  to  achieve  this 
validation,  image  reconstruction  was  performed  for  the  full-body  in-vivo  mouse  OAT  imaging  in 
addition  of  the  numerical  phantom  study.  The  results  were  reported  in  an  oral  presentation  for 
SPIE  2015.  The  schematic  of  the  imaging  module  is  shown  in  fig.  2.  The  comparison  of  mouse 
images  using  unmatched  imaging  model  to  the  matched  imaging  model  is  shown  in  fig.  2.  It  can 
be  seen  clearly  that  using  the  unmatched  backprojection  imaging  model  has  not  compromised  the 
resolution  of  the  reconstructed  images.  This  validation  is  pivotal  in  achieving  the  goal  of  iterative 
image  reconstruction  of  breast  imaging.  For  more  detail  please  refer  to  Appendix  IV. 


3b.  (i)  GPU  implementation  of  image  reconstruction  algorithms:  Improved  GPU-based 
implementations  of  a  numerical  imaging  model  and  its  adjoint  have  been  developed  for  use  with 
general  gradient-based  iterative  image  reconstruction  algorithms.  Particularly,  two  types  of 
computation-reduced  discretization  methods  have  been  employed;  a  parallel  fast  GPU-based 
Fourier  transform  (FFT)  algorithm  was  employed  to  accelerate  the  calculation  of  the  temporal 
convolution  with  ultrasonic  transducer  responses;  and  a  volume-reduction  method  is  proposed  to 
reduce  the  computation  for  applications  with  irregular  field-of-view  (breast  imaging).  The  results 
suggest  that  the  proposed  implementation  is  more  than  five  times  faster  than  previous 
implementations  for  a  single  GPU.  In  addition,  the  algorithm  has  also  been  developed  to  use 
multiple  GPUs  further  reducing  the  computational  time.  The  work  was  presented  in  SPIE 
Photonics  West,  2015. 


(ii)  Development  of  Innovative  Forward  Projection  numerical  scheme  to  accomplish  Breast 
Imaging:  The  unmatched  back-projection  imaging  model,  and  GPU  implementation  are  two 
major  accomplishments  during  the  first  year  of 
developing  imaging  model  and  optimization 
algorithm.  The  next  task  is  to  be  able  to  perform 
image  reconstruction  for  the  breast.  The  iterative 
image  reconstruction  involved  as  much  as  up  to 
35000  integration  step  over  a  volume  of  up  to  70 
mm  x  70  mm  x  90  mm,  with  the  discretization 
requirement  much  less  than  1  mm.  The  iterative 
reconstruction  tomography  required  more 
innovation.  To  achieve  even  faster  images,  I 
exploited  the  symmetry  of  the  problem.  The 
schematic  of  the  object  volume  is  shown  in  Fig. 

In  previous  implementation  the  object  volume 
was  fully  contained  inside  a  sphere.  The 
integration  was  performed  over  the  surface  of  Figure  3:  Schematic  of  the  integration  volume  for  the  forward 
intersection  between  the  two  spheres  (one  of  the  projection  of  imaging  model. 

object  and  second  defined  by  the  transducers’  location  to  receive  pressure  wave  at  a  given  time) 
to  calculate  the  forward  projection  (blue  and  grey  shaded  regions).  However,  this  was  infeasible 
for  breast  imaging  application.  In  this  case,  I  devised  the  hemispherical  object  volume  (red 
shaded  region)  and,  the  chest  wall  was  modeled  by  a  plane.  I  derived  analytically  the  intersection 
between  the  plane  and  the  integration  surface  to  eliminate  undesirable  integration  steps.  This 


Surface  used  in  integral 
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allowed  us  to  conduct  studies  for  numerical  breast  phantoms  of  realistic  sizes  and  clinical 
studies. 

Task  4:  Validate  prototype  imager  and  image  reconstruction  algorithms 

4a.  Experimental  phantom  studies: 

We  conducted  experimental  phantom  studies  using  the  first- 
generation  prototype  OAT  breast  imaging  system  built  by  our 
collaborators  at  TomoWave  Inc.  The  first  experiment  uses  a 
gelatin  phantom  mimicking  breast  tissue.  A  gelatin  phantom 
was  created  to  mimic  normal  breast  tissue,  with  an  embedded 
gelatin  sphere  to  represent  tumor  and  two  crossing  tubes  filled 
with  an  optical-absorbing  fluid  to  represent  blood  vessel.  The 
phantom  was  imaged.  The  pressure  data  was  collected  for 
1536  time  samples  with  a  sampling  frequency  of  12.5  MHz  at 
1800  tomographic  views.  The  reconstructed  OAT  images  are 
displayed  in  Fig.  4.  These  images  accurately  revealed  the  phantom  structures. 

4b.  Use  of  clinical  image  data: 

We  also  conducted  preliminary  that  employed  clinical 
data.  Figure  5  shows  an  OAT  image.  Pressure  data  is 
collected  at  only  300  tomographic  views.  It  is  very 
encouraging  to  observe  that  some  major  vascular 
structures  can  be  clearly  reconstructed  along  with  some 
peripheral  small  capillary  vessels.  The  study  helped 
immensely  for  phase-II  design  of  imaging  module  and 
imaging  reconstruction  methods. 

4.  KEY  RESEARCH  ACCOMPLISHMENTS 


Figure  5:  Maximum  intensity  projection 
image  of  a  female  breast  (clinical  study). 


Figure  4:  Maximum  intensity 
projection  image  of  an  experimental 
phantom  contains  tubes  and  a  sphere. 


•  A  novel  OAT  imaging  system  for  breast  imaging  was  designed  and  implemented. 

•  Image  reconstruction  methods  were  developed,  implemented,  and  tested. 

•  The  developed  system  and  algorithms  were  employed  in  a  pilot  clinical  study  in  which 
3D  OAT  images,  in  vivo,  were  acquired.  The  images  revealed  the  presence  of  vessels 
and  other  blood  filled  structures. 

5.  Conclusion 

This  project  involved  the  development  and  clinical  validation  of  an  innovative  3D  breast  imaging 
system  that  provides  co-registered  functional  and  anatomical  information  without  the  use  of 
ionizing  radiation.  The  system  is  first  of  its  kind  to  produce  comprehensive  clinical  information 
based  on  co-registered  volumetric  images  of  two  types  of  OAT  images  depicting  the  absorbed 
optical  energy  (blood  distribution  and  its  oxygenation)  and  three  types  of  UST  images  (speed-of- 
sound,  ultrasonic  reflectance,  and  ultrasonic  attenuation  distributions)  of  the  breast  tissue.  This 
will  provide  an  extremely  rich  set  of  complementary  anatomical  and  functional  3D  information 
to  the  radiologist,  enabling  early  detection  of  small  tumors  (such  as  in  situ  lesions),  accurate 
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staging  of  more  advanced  tumors,  and  diagnostic  characterization  of  the  tumor  potential  for 
aggressive  growth  and  metastatic  activity.  Early  tumors  (<  6  mm)  are  undetectable  with  the 
existing  screening  mammography,  especially  in  a  dense  breast  of  younger  women.  Unlike  X-ray 
based  or  nuclear  imaging  methods,  the  proposed  imaging  method  is  radiation-safe,  pain- free  (no 
compression  required),  toxicity-free,  and  much  less  expensive.  The  new  ultrasound  plus 
optoacoustic  tomography  platform  is  expected  to  transform  breast  cancer  care  by  providing  a 
better  screening  method  and  serve  as  a  highly  valuable  tool  for  monitoring  and  identifying 
effective  cancer  therapies. 
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Authors:  Fatima  Anis,  Yang  Lou,  Andre  Conjusteau,  Sergey  Ermilov,  Alexander  Oraevsky  and 
Mark  A.  Anastasio 

Conference:  SPIE  Photonics  West-  2014;8943,  San  Francisco  CA 

2b-  Title:  Investigation  of  a  method  for  laser-induced  ultrasound  tomography  that  eliminates  the 
need  for  ray  tracing 

Authors:  Fatima  Anis,  Yang  Lou,  Andre  Conjusteau,  Sergey  Ermilov, 

Alexander  Oraevsky  and  Mark  A.  Anastasio 

Conference:  SPIE  Medical  Imaging  2015:  Ultrasonic  Imaging  and  Tomography 

3b-  Title:  Accelerated  iterative  image  reconstruction  in  three-dimensional  optoacoustic 
tomography 

Authors:  Fatima  Anis,  Yang  Lou,  Kun  Wang,  Richard  Su,  Tanmayi  Oruganti,  Andre 
Conjusteau,  Sergey  Ermilov,  Alexander  A.  Oraevsky  and  Mark  A.  Anastasio 
Conference:  SPIE  Photonics  West,  2015 

4b-  Title:  Waveform  Inversion  with  Source  Encoding  for  Breast  Speed-of-Sound  Reconstruction 
in  Ultrasound  Computed  Tomography 

Authors:  Kun  Wang,  Thomas  Mathew,  Fatima  Anis,  Cuiping  Li,  Neb  Duric,  and  Mark  A. 
Anastasio 

Conference:  SPIE  Medical  Imaging  2015:  Ultrasonic  Imaging  and  Tomography 
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5b-  Title:  Three-dimensional  laser  optoacoustic  and  laser  ultrasound  imaging  system  for 
biomedical  research 

Authors:  Sergey  A.  Ermilov,  Rishard  Su,  Andre  Conjusteau,  Kun  Wag,  Fatima  Anis,  Mark  A. 
Anastasio,  Alexander  A.  Oreavsky 

Conference:  SPIE:  Photon  Plus  Ultrasound:  Imaging  and  Sensing  2015;  9323 


7.  Training 

The  year  has  been  very  fruitful  as  I  continue  to  benefit  from  many  scholarly  activities  and  kept 
adding  to  my  skills  as  imaging  scientist. 

1-  Attending  course  (E62  BME  500  67):  In  Fall  2014, 1  participated  in  BME  course  on  the 
imaging  science.  The  class  held  weekly  for  two  hours.  Some  of  the  study  topics  included 
continuous  and  discrete  object  representations,  imaging  operators,  image  statistics, 
imaging  quality  assessment,  ideal  observer,  Hotelling  observers  and  imaging  errors.  The 
course  provided  me  with  a  formal  training  as  biomedical  image  scientist. 

2-  Conferences: 

i-  SPIE  Photonics  West,  2014  San  Francisco  CA:  The  conference  is  the  major 
international  meeting  held  annually  for  the  optoacoustic/photoacpustic  imaging.  I 
attended  this  meeting  and  was  greatly  benefited  from  the  presentation  and  poster  sessions 
as  well  as  constructive  meetings  with  the  Prof.  Oraevsky  and  his  team  about  the  hybrid 
OAT/  USCT  breast  imager. 

ii-  SPIE  Medical  Imaging,  2014,  San  Diago  CA:  This  conference  is  another  major 
international  conference,  which  holds  annually  and  covers  broad  range  of  topics 
concerning  medical  imaging  and  diagnostics.  I  was  specially  benefited  from  many  talks 
and  poster  presentations  about  the  image  quality  assessment.  Moreover,  the  dedicated 
Ultrasonic  Imaging  and  Tomography  sessions  on  ultrasound  imaging  provided  me  with  a 
great  opportunity  to  learn  about  ultrasound  tomography  in  medical  imaging. 

iii-  SPIE  Photonics  West,  2015  San  Francisco  CA 

3-  Visiting  TomoWave  Inc,  Houston  TX: 

I  visited  TomoWave  Inc.  in  December  2013.  During  the  visit,  I  benefited  from  lab  tours 
and  learned  about  practical  aspects  of  OAT/USCT  imaging. 

4-  Algorithm  development  and  GPU  computing:  I  continued  to  establish  more  skills 
towards  algorithm  development.  One  of  the  major  achievements  towards  this  end  was  to 
learn  CUDA  programming  from  other  group  members.  This  training  will  continue  to 
benefit  me  through  the  remainder  of  the  project. 

5-  Image  rendering  software:  I  learned  image  rendering  softwares  for  3D  visualization  of 
the  imaging.  This  is  important  and  allows  one  to  draw  inference  form  the  reconstructed 
images. 
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8.  Reportable  Outcomes 

Nothing  to  report 

9.  Other  Achievements 

Nothing  to  report 
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Waveform  Inversion  With  Source  Encoding 
for  Breast  Sound  Speed  Reconstruction 
in  Ultrasound  Computed  Tomography 

Kun  Wang,  Member,  IEEE ,  Thomas  Matthews,  Student  Member,  IEEE ,  Fatima  Anis,  Cuiping  Li, 
Neb  Duric,  and  Mark  A.  Anastasio,  Senior  Member,  IEEE 


Abstract — Ultrasound  computed  tomography  (USCT)  holds 
great  promise  for  improving  the  detection  and  management 
of  breast  cancer.  Because  they  are  based  on  the  acoustic  wave 
equation,  waveform  inversion-based  reconstruction  methods 
can  produce  images  that  possess  improved  spatial  resolution 
properties  over  those  produced  by  ray-based  methods.  How¬ 
ever,  waveform  inversion  methods  are  computationally  de¬ 
manding  and  have  not  been  applied  widely  in  USCT  breast 
imaging.  In  this  work,  source  encoding  concepts  are  employed 
to  develop  an  accelerated  USCT  reconstruction  method  that 
circumvents  the  large  computational  burden  of  conventional 
waveform  inversion  methods.  This  method,  referred  to  as  the 
waveform  inversion  with  source  encoding  (WISE)  method,  en¬ 
codes  the  measurement  data  using  a  random  encoding  vector 
and  determines  an  estimate  of  the  sound  speed  distribution  by 
solving  a  stochastic  optimization  problem  by  use  of  a  stochas¬ 
tic  gradient  descent  algorithm.  Both  computer  simulation  and 
experimental  phantom  studies  are  conducted  to  demonstrate 
the  use  of  the  WISE  method.  The  results  suggest  that  the 
WISE  method  maintains  the  high  spatial  resolution  of  wave¬ 
form  inversion  methods  while  significantly  reducing  the  com¬ 
putational  burden. 


I.  Introduction 

After  decades  of  research  [l]-[4],  advancements  in 
hardware  and  computing  technologies  are  now  fa¬ 
cilitating  the  clinical  translation  of  ultrasound  computed 
tomography  (USCT)  for  breast  imaging  applications  [2], 
[5]— [8] .  USCT  holds  great  potential  for  improving  the 
detection  and  management  of  breast  cancer  because  it 
provides  novel  acoustic  tissue  contrasts,  is  radiation-  and 
breast-compression- free,  and  is  relatively  inexpensive  [9], 
[10].  Several  studies  have  reported  the  feasibility  of  USCT 
for  characterizing  breast  tissues  [2],  [4]— [6] ,  [10],  [11].  Al¬ 
though  some  USCT  systems  are  capable  of  generating 
three  images  that  depict  the  breast’s  acoustic  reflectivity, 


Manuscript  received  October  20,  2014;  accepted  December  28,  2014. 
This  work  was  supported  in  part  by  NIH  awards  EB010049,  CA1744601, 
EB01696301,  and  DOD  Award  US  ARMY  W81XWH-13-1-0233. 

K.  Wang,  T.  P.  Matthews,  F.  Anis,  and  M.  A.  Anastasio  are  with  the 
Department  of  Biomedical  Engineering,  Washington  University  in  St. 
Louis,  St.  Louis,  MO  63130,  USA  (e-mail:  anastasio@wustl.edu). 

C.  Li  and  N.  Duric  are  with  Delphinus  Medical  Technologies,  Plym¬ 
outh,  MI  48170,  USA. 

N.  Duric  is  also  with  Karmanos  Cancer  Institute,  Wayne  State  Univer¬ 
sity,  Detroit,  MI  48201,  USA. 

DOI  http://dx.doi.org/10.1109/TUFFC.2014.006788 


acoustic  attenuation,  and  sound  speed  distributions,  this 
study  will  focus  on  the  reconstruction  of  the  sound  speed 
distribution. 

A  variety  of  USCT  imaging  systems  have  been  devel¬ 
oped  for  breast  sound  speed  imaging  [5],  [7],  [10],  [12]— [15] . 
In  a  typical  USCT  experiment,  acoustic  pulses  that  are 
generated  by  different  transducers  are  employed,  in  turn, 
to  insonify  the  breast.  The  resulting  wavefield  data  are 
measured  by  an  array  of  ultrasonic  transducers  that  are 
located  outside  of  the  breast.  Here  and  throughout  the 
manuscript,  a  transducer  that  produces  an  acoustic  pulse 
will  be  referred  to  as  an  emitter;  the  transducers  that 
receive  the  resulting  wavefield  data  will  be  referred  to  as 
receivers.  From  the  collection  of  recorded  wavefield  data, 
an  image  reconstruction  method  is  utilized  to  estimate  the 
sound  speed  distribution  within  the  breast  [5],  [7],  [10]. 

The  majority  of  USCT  image  reconstruction  methods 
for  breast  imaging  investigated  to  date  have  been  based 
on  approximations  to  the  acoustic  wave  equation  [12], 
[16]— [24] .  A  relatively  popular  class  of  methods  is  based 
on  geometrical  acoustics.  Such  methods  are  commonly  re¬ 
ferred  to  as  ray-based  methods.  These  methods  involve 
two  steps.  First,  time-of- flight  (TOF)  data  corresponding 
to  each  emitter-receiver  pair  are  estimated  [25].  Under  a 
geometrical  acoustics  approximation,  the  TOF  data  are 
related  to  the  sound  speed  distribution  via  an  integral  ge¬ 
ometry,  or  ray-based,  imaging  model  [16],  [26].  Second,  by 
use  of  the  measured  TOF  data  and  the  ray-based  imaging 
model,  a  reconstruction  algorithm  is  employed  to  estimate 
the  sound  speed  distribution.  Although  ray-based  meth¬ 
ods  can  be  computationally  efficient,  the  spatial  resolu¬ 
tion  of  the  images  they  produce  is  limited  due  to  the  fact 
that  diffraction  effects  are  not  modeled  [23],  [27].  This  is 
undesirable  for  breast  imaging  applications,  in  which  the 
ability  to  resolve  fine  features  (e.g.,  tumor  spiculations)  is 
important  for  distinguishing  healthy  from  diseased  tissues. 

USCT  reconstruction  methods  based  on  the  acoustic 
wave  equation,  also  known  as  full-wave  inverse  scattering 
or  waveform  inversion  methods,  have  also  been  explored 
for  a  variety  of  applications  including  medical  imaging 
[12],  [22],  [23],  [28]  and  geophysics  [29]— [31] .  Because  they 
account  for  higher-order  diffraction  effects,  waveform  in¬ 
version  methods  can  produce  images  that  possess  high¬ 
er  spatial  resolution  than  those  produced  by  ray-based 
methods  [23],  [28].  However,  conventional  waveform  inver¬ 
sion  methods  are  iterative  in  nature  and  require  the  wave 
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equation  to  be  solved  numerically  a  large  number  of  times 
at  each  iteration.  Consequently,  such  methods  can  be  ex¬ 
tremely  computationally  burdensome.  For  special  geom¬ 
etries  [12],  [32],  efficient  numerical  wave  equation  solvers 
have  been  reported.  However,  apart  from  special  cases,  the 
large  computational  burden  of  waveform  inversion  meth¬ 
ods  has  hindered  their  widespread  application. 

A  natural  way  to  reduce  the  computational  complexity 
of  the  reconstruction  problem  is  to  reformulate  it  in  a  way 
that  permits  a  reduction  in  the  number  of  times  the  wave 
equation  needs  to  be  solved.  In  the  geophysics  literature, 
source  encoding  methods  have  been  proposed  to  achieve 
this  [29]— [31] .  When  source  encoding  is  employed,  at  each 
iteration  of  a  prescribed  reconstruction  algorithm,  all  of 
the  acoustic  pulses  produced  by  the  emitters  are  combined 
(or  encoded)  by  use  of  a  random  encoding  vector.  The 
measured  wavefield  data  are  combined  in  the  same  way. 
As  a  result,  the  wave  equation  may  need  to  be  solved  as 
few  as  twice  at  each  algorithm  iteration.  In  conventional 
waveform  inversion  methods,  this  number  would  be  equal 
to  twice  the  number  of  emitters  employed.  Although  con¬ 
ventional  waveform  inversion  methods  may  require  fewer 
algorithm  iterations  to  obtain  a  specified  image  accuracy 
compared  with  source  encoded  methods,  as  demonstrated 
later,  the  latter  can  greatly  reduce  the  overall  number  of 
times  the  wave  equation  needs  to  be  solved. 

In  this  study,  a  waveform  inversion  with  source  encod¬ 
ing  (WISE)  method  for  USCT  sound  speed  reconstruction 
is  developed  and  investigated  for  breast  imaging  with  a 
circular  transducer  array.  The  WISE  method  determines 
an  estimate  of  the  sound  speed  distribution  by  solving 
a  stochastic  optimization  problem  by  use  of  a  stochastic 
gradient  descent  algorithm  [30],  [33].  Unlike  previously 
studied  waveform  inversion  methods  that  were  based  on 
the  Helmholtz  equation  [22] ,  [23] ,  the  WISE  method  is  for¬ 
mulated  by  use  of  the  time-domain  acoustic  wave  equation 
[34] -[36]  and  uses  broad-band  measurements.  The  wave 
equation  is  solved  by  use  of  a  computationally  efficient 
k-space  method  that  is  accelerated  by  use  of  graphics  pro¬ 
cessing  units  (GPUs).  To  mitigate  the  interference  of  the 
emitter  on  its  neighboring  receivers,  a  heuristic  data  re¬ 
placement  strategy  is  proposed.  The  method  is  validated 
in  computer  simulation  studies  that  include  modeling  er¬ 
rors  and  other  physical  factors.  The  practical  applicability 
of  the  method  is  further  demonstrated  in  studies  involving 
experimental  breast  phantom  data. 

The  remainder  of  the  paper  is  organized  as  follows.  In 
Section  II,  USCT  imaging  models  in  their  continuous  and 
discrete  forms  are  reviewed.  A  conventional  waveform  in¬ 
version  method  and  the  WISE  method  for  sound  speed  re¬ 
construction  are  formulated  in  Section  III.  The  computer 
simulation  studies  and  corresponding  numerical  results  are 
presented  in  Sections  IV  and  V,  respectively.  In  Section 
VI,  the  WISE  method  is  further  validated  in  experimental 
breast  phantom  studies.  Finally,  the  paper  concludes  with 
a  discussion  in  Section  VII. 


II.  Background:  USCT  Imaging  Models 

In  this  section,  imaging  models  that  provide  the  ba¬ 
sis  for  image  reconstruction  in  waveform  inversion-based 
USCT  are  reviewed. 

A.  USCT  Imaging  Model  in  Its  Continuous  Form 

Although  a  digital  imaging  system  is  properly  described 
as  a  continuous-to-discrete  (C-D)  mapping  (chapter  7  in 
[37]),  for  simplicity,  a  USCT  imaging  system  is  initially 
described  in  its  continuous  form  below. 

In  USCT  breast  imaging,  a  sequence  of  acoustic  pulses 
is  transmitted  through  the  breast.  We  denote  each  acous¬ 
tic  pulse  by  sm( r,t)  E  L2(M3  x  [0,oo)),  where  each  pulse  is 
indexed  by  an  integer  m  for  m  =  0,  1,  ...,  M  —  1  with  M 
denoting  the  total  number  of  acoustic  pulses.  Although  it 
is  spatially  localized  at  the  emitter  location,  each  acoustic 
pulse  can  be  expressed  as  a  function  of  space  and  time. 
When  the  rath  pulse  propagates  through  the  breast,  it 
generates  a  pressure  wavefield  distribution  denoted  by 
pm( r,£)  E  L2(IR3  x  [0,oo)).  If  acoustic  absorption  and  mass 
density  variations  are  negligible,  pm(r,t)  in  an  unbounded 
medium  satisfies  the  acoustic  wave  equation  [38]: 

1  d 2 

V2pm( T,t)  -  2,  ,  a.2PmM  =  -4ns Jr, t),  (1) 

c  ( r)  at 

where  c(r)  is  the  sought-after  sound  speed  distribution. 
Eq.  (1)  can  be  expressed  in  operator  form  as 

Pm(r,t)  =  Hcsjr,t),  (2) 

where  the  linear  operator  Ti c  :  L2(IR3  x  [0,oo))  I-V 
L2(R3  x[0,oo))  denotes  the  action  of  the  wave  equation 
and  is  independent  of  the  index  of  ra.  The  superscript  c 
indicates  the  dependence  of  Tic  on  c(r). 

Consider  that  pm(r,t)  is  recorded  outside  of  the  object 
for  r  E  and  t  E  [0,7],  where  Qm  C  M3  denotes  a  con¬ 
tinuous  measurement  aperture.  In  this  case,  when  discrete 
sampling  effects  are  neglected,  the  imaging  model  can  be 
described  as  a  continuous-to-continuous  mapping  as 

=  MmHcsJr,t),  for  m  =  0,1  ,  —  ,M  -  1,  (3) 

where  gjr,t)  €  L2(T2m  x  [0,T])  denotes  the  measured  data 
function  and  the  operator  J\4m  is  the  restriction  of  Hc  to 
x  [0,  T].  The  ra-dependent  operator  A im  allows  (3)  to 
describe  USCT  imaging  systems  in  which  the  measure¬ 
ment  aperture  varies  with  emitter  location.  Here  and 
throughout  the  manuscript,  we  will  refer  to  the  process  of 
firing  one  acoustic  pulse  and  acquiring  the  corresponding 
wavefield  data  as  one  data  acquisition  indexed  by  ra.  The 
USCT  reconstruction  problem  in  its  continuous  form  is  to 
estimate  the  sound  speed  distribution  c(r)  by  use  of  (3) 
and  the  data  functions  {gm(r,t)}^lQ. 
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B.  USCT  Imaging  Model  in  Its  Discrete  Forms 

A  digital  imaging  system  is  accurately  described  by  a 
continuous-to-discrete  (C-D)  imaging  model,  which  is  typ¬ 
ically  approximated  in  practice  by  a  discrete-to-discrete 
(D-D)  imaging  model  to  facilitate  the  application  of  itera¬ 
tive  image  reconstruction  algorithms.  A  C-D  description 
of  the  USCT  imaging  system  is  provided  in  Appendix  A. 
Below,  a  D-D  imaging  model  for  waveform-based  USCT 
is  presented.  This  imaging  model  will  be  employed  subse¬ 
quently  in  the  development  of  the  WISE  method  in  Sec¬ 
tion  III. 

Construction  of  a  D-D  imaging  model  requires  the  in¬ 
troduction  of  finite-dimensional  approximate  representa¬ 
tions  of  the  functions  c(r)  and  sm(r,t),  which  will  be  de¬ 
noted  by  the  vectors  c  G  and  sm  G  RNL.  Here,  N  and 
L  denote  the  number  of  spatial  and  temporal  samples, 
respectively,  employed  by  the  numerical  wave  equation 
solver.  In  waveform-based  USCT,  the  way  in  which  c(r) 
and  sm(r,t)  are  discretized  to  form  c  and  sm  is  dictated  by 
the  numerical  method  employed  to  solve  the  acoustic  wave 
equation.  In  this  study,  we  employ  a  pseudospectral  k- 
space  method  [34]  [36].  Accordingly,  c(r)  and  sm(r,t)  are 
sampled  on  Cartesian  grid  points  as 


characteristics,  such  as  finite  aperture  size  and  temporal 
delays.  For  simplicity,  we  assume  that  the  transducers  are 
point-like  in  this  study.  When  the  receiver  and  grid  point 
locations  do  not  coincide,  an  interpolation  method  is  re¬ 
quired.  As  an  example,  when  a  nearest-neighbor  interpola¬ 
tion  method  is  employed,  the  elements  of  Mm  are  defined 
as 


\M-m\n™L+l,nL+l 


jl,  fc>r  n  =  ^m(^rec)5 

[0,  otherwise, 


(7) 


where  [^-mln^L+inL+i  denotes  the  element  of  Mm  at  the 
(nrecL  +  Z)th  row  and  the  (nL  +  Z)th  column,  and  Tm(nTec) 
denotes  the  index  of  the  grid  point  that  is  closest  to 
r(m,nrec).  Here,  r(m,7irec)  denotes  the  location  of  the  nrecth 
receiver  in  the  mth  data  acquisition.  Eq.  (6)  represents  the 
D-D  imaging  model  that  will  be  employed  in  the  remain¬ 
der  of  this  study. 


III.  Waveform  Inversion  With  Source 
Encoding  for  USCT 

A.  Sequential  Waveform  Inversion  in  Its  Discrete  Form 


[c]n  =  c(r„),  and  [sm]nL+l  =  sm{rnlW), 

.  n  =  0,l,-,N-l  (4) 

tor  , 

l  =  0,1, -,L  — 1 

where  Af  denotes  the  temporal  sampling  interval  and  rn 
denotes  the  location  of  the  nth  point. 

For  a  given  c  and  sm,  the  pseudospectral  k-space  meth¬ 
od  can  be  described  in  operator  form  as 

P™  =  Hcsm,  (5) 

where  the  matrix  Hc  is  of  dimension  NL  x  NL  and  repre¬ 
sents  a  discrete  approximation  of  the  wave  operator  Hc 
defined  in  (2),  and  the  vector  represents  the  estimated 
pressure  data  at  the  grid  point  locations  and  has  the  same 
dimension  as  sm.  The  superscript  a  indicates  that  these 
values  are  approximate,  i.e.,  [p^lnL+i  ~  pm(rn,ZA*).  We  re¬ 
fer  the  readers  to  [34] -[36]  for  additional  details  regarding 
the  pseudospectral  k-space  method. 

Because  the  pseudospectral  k-space  method  yields  sam¬ 
pled  values  of  the  pressure  data  on  a  Cartesian  grid,  a 
sampling  matrix  Mm  is  introduced  to  model  the  USCT 
data  acquisition  process  as 

gS.  =  =  MmHcsm.  (6) 

Here,  the  NrecL  x  NL  sampling  matrix  Mm  extracts  the 
pressure  data  corresponding  to  the  receiver  locations  on 
the  measurement  aperture  Um,  with  Nrec  denoting  the 
number  of  receivers.  The  vector  denotes  the  predicted 
data  that  approximates  the  true  measurements.  In  prin¬ 
ciple,  Mm  can  be  constructed  to  incorporate  transducer 


A  conventional  waveform  inversion  method  that  does 
not  utilize  source  encoding  will  be  employed  as  a  reference 
for  the  developed  WISE  method  and  is  briefly  described 
below.  Like  other  conventional  approaches,  this  method 
sequentially  processes  the  data  acquisitions  {gm}  (in  an 
arbitrary  order)  for  m  —  0,  1,  ...,  M  —  1  at  each  iteration 
of  the  associated  algorithm.  As  such,  we  will  refer  to  the 
conventional  method  as  a  sequential  waveform  inversion 
method. 

A  sequential  waveform  inversion  method  can  be  formu¬ 
lated  as  a  nonlinear  numerical  optimization  problem: 

c  =  arg  min{jF(c)  +  /371(c)},  (8) 

c 

where  F(c),  7£(c),  and  (3  denote  the  data  fidelity  term,  the 
penalty  term,  and  the  regularization  parameter,  respec¬ 
tively.  The  data  fidelity  term  jF(c)  is  defined  as  a  sum  of 
squared  1 2  norms  of  the  data  residuals  corresponding  to 
all  data  acquisitions  as 

M— 1 

T(c)  =  2E  II  8™  -  MmHcsm  II2,  (9) 

771  =  0 

where  gm  G  WN  e°L  denotes  the  measured  data  vector  at 
the  mth  data  acquisition.  The  choice  of  the  penalty  term 
will  be  addressed  in  Section  IV. 

The  gradient  of  F(c)  with  respect  to  c,  denoted  by  J, 
will  be  computed  by  discretizing  an  expression  for  the 
Frechet  derivative  that  is  derived  assuming  a  continuous 
form  of  (9).  The  Frechet  derivative  is  described  in  Appen¬ 
dix  B.  Namely,  the  gradient  is  approximated  as 
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M—l 

[J]n  =  y^[^m]n 

m= 0 

M  — 17-2 

~  m3  ^Yy^m\nL+(L-l)  (10) 

Tin  m=0l=l 

^  \P  m\nL+l  —  1  ~  ^[Pm]n7  +  /  [Pm]n7  +  /  +  1 

A1 

where  Jm  denotes  the  gradient  of  ||  gm  —  MmHcsm  ||2/2 
with  respect  to  c  and  the  vector  contains  samples  that 
approximate  adjoint  wavefield  qm(r,t)  that  satisfies  (34)  in 
Appendix  B.  By  use  of  the  psendospectral  k-space  meth¬ 
od,  can  be  calculated  as 


tion  solver.  Note  that  it  appears  in  Lines  6,  7,  and  11. 
Because  Lines  6  and  7  have  to  be  executed  M  times  to 
process  all  of  the  data  acquisitions  and  the  line  search 
method  must  pass  through  all  M  emitters  at  least  once, 
the  wave  equation  solver  has  to  be  executed  at  least  3 M 
times  at  each  algorithm  iteration.  The  line  search  in  Line 
11  searches  for  a  step  size  along  the  direction  of  —  J  so 
that  the  cost  function  is  reduced  by  use  of  a  classic  trial- 
and-error  approach  [39].  Note  that,  in  general,  the  line 
search  will  require  more  than  M  applications  of  Hc,  so 
3 M  represents  a  lower  bound  on  the  total  number  of  wave 
equation  solver  runs  per  iteration. 

B.  Stochastic  Optimization- Based  WISE 


9m  =  (11) 

where 

fem  “  ^m\x-\n)L+{L-iy  ^ n  ^  ^2) 

0,  otherwise. 

Here,  Nm  =  {n  :  Xm(nrec),nvec  =  0,1,-, Nrec  -  1},  and 
X~l  denotes  the  inverse  mapping  of  Xm. 

Given  the  explicit  form  of  J  in  (10),  a  variety  of  opti¬ 
mization  algorithms  can  be  employed  to  solve  (8)  [39]. 
Algorithm  1  describes  a  gradient  descent-based  sequential 
waveform  inversion  method.  Note  that  at  every  algorith¬ 
mic  iteration,  the  sequential  waveform  inversion  method 
updates  the  sound  speed  estimate  only  once  using  the  gra¬ 
dient  J  accumulated  over  all  Jm  for  m  =  0,  1,  ...,  M  —  1. 
This  is  unlike  the  Kaczmarz  method — also  known  as  the 
algebraic  reconstruction  technique  [16],  [19],  [40] — that 
updates  the  sound  speed  estimate  multiple  times  in  one 
algorithmic  iteration.  In  Line  10  of  Algorithm  1,  JR  de¬ 
notes  the  gradient  of  E(c)  with  respect  to  c. 

Algorithm  1:  Gradient  descent-based  sequential  wave¬ 
form  inversion. 

Input:  {gm},{sm},c(0) 

Output:  c 

•  1:  k  0  {k  is  the  number  of  algorithm  iteration.} 

•  2:  while  stopping  criterion  is  not  satisfied  do 

•  3:  k  <—  k  +  1 

•  4:  J<-0 

•  5:  for  m  :=  0  to  M  —  1  do 

•  6:  <-  Hcsm  {m  is  the  index  of  the  emitter.} 

•  7:  q^  Hcrm  {rm  is  calculated  via  (12).} 

•  8:  J  *  J  T  J m  Urn  is  calculated  via  (10).} 

•  9:  end  for 

•  10:  J  <—  J  +  /3Jr 

•  11:  Determine  step  size  A  via  a  line  search 

•  12:  c®  <-  c^-1)  -  AJ 

•  13:  end  while 

•  14:  c  =  cw 

In  Algorithm  1,  Hc  is  the  most  computationally  bur¬ 
densome  operator,  representing  one  run  of  the  wave  equa¬ 


To  alleviate  the  large  computational  burden  presented 
by  sequential  waveform  inversion  methods  (e.g.,  Algorithm 
1),  a  source  encoding  method  has  been  proposed  [22],  [29], 
[41].  This  method  has  been  formulated  as  a  stochastic  op¬ 
timization  problem  and  solved  by  various  stochastic  gra¬ 
dient-based  algorithms  [30],  [31].  In  this  section,  we  adapt 
the  stochastic  optimization-based  formulation  in  [30]  to 
find  a  solution  of  (8). 

Algorithm  2.  WISE  Method. 

Input:  {g4,{sm},c(0) 

Output:  c 


•  1:  k  <—  0  {k  is  the  number  of  algorithm  iteration.} 

•  2:  while  stopping  criterion  is  not  satisfied  do 

•  3:  k  <—  k  +  1 


•  4:  Draw  elements  of  w  from  independent  and  identi¬ 
cal  Rademacher  distribution. 


pw  Hcsw  {sw  is  calculated  via  (14).} 
qw  <—  H crw  {rw  is  calculated  via  (17).} 

J  <—  Jw  +  (3Jr  {Jw  is  calculated  via  (16).} 

Determine  step  size  A  by  use  of  line  search 
c(k)  c(k- 1)  _  yj 

•  10:  end  while 

•  11:  c  = 


•  5 

•  6 

•  7 

•  8 
•  9 


The  WISE  method  seeks  to  minimize  the  same  cost 
function  as  the  sequential  waveform  inversion  method, 
namely  (8).  However,  to  accomplish  this,  the  data  fidelity 
term  in  (9)  is  reformulated  as  the  expectation  of  a  random 
quantity  as  [29]— [31] ,  [33],  [41],  [42] 

U s(c)  =  Ew{i  ||  gw  —  MHcsw  ||2},  (13) 

where  Ew  denotes  the  expectation  operator  with  respect 
to  the  random  source  encoding  vector  w  £  IRM,  M  =  Mm 
is  the  sampling  matrix  that  is  assumed  to  be  identical  for 
m  —  0,  1,  ...,  M  —  1,  and  gw  and  sw  denote  the  w-encoded 
data  and  source  vectors,  defined  as 


M—l 


M—l 


T  =  and  sW  =  Etw]>»s>»’  (14) 


m= 0 


m= 0 
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Fig.  1.  Schematic  of  a  USCT  system  with  a  circular  transducer  array 
whose  elements  are  indexed  from  0  to  255.  It  shows  the  first  data  acquisi¬ 
tion,  where  element  0  (in  red)  is  emitting  an  acoustic  pulse,  whereas  all 
256  elements  are  receiving  signals.  The  region  of  interest  (ROI)  is  shaded 
in  gray,  and  the  dashed  square  box  represents  the  physical  dimensions 
(128  x  128  mm2)  of  all  reconstructed  images. 


[Tw]nL+l  =  hMpW  £  b-»£+(£-0’  lfnGN>  .  (17) 

[0,  otherwise 

Here,  we  drop  the  subscript  m  of  both  T~\n)  and  N  be¬ 
cause  we  assume  M  to  be  identical  for  all  data  acquisi¬ 
tions.  Various  probability  density  functions  have  been  pro¬ 
posed  to  describe  the  source  encoding  vector  w  [29],  [31], 
[41].  In  this  study,  we  employed  a  Rademacher  distribu¬ 
tion  as  suggested  by  [29] ,  in  which  case  each  element  of  w 
had  a  50%  chance  of  being  either  +1  or  —1. 


IV.  Description  of  Computer  Simulation  Studies 

Two-dimensional  computer  simulation  studies  were 
conducted  to  validate  the  WISE  method  for  breast  sound 
speed  imaging  and  demonstrate  its  computational  ad¬ 
vantage  over  the  standard  sequential  waveform  inversion 
method. 


respectively.  It  has  been  demonstrated  that  (9)  and  (13) 
are  mathematically  equivalent  when  w  possesses  a  zero 
mean  and  an  identity  covariance  matrix  [30],  [33],  [42].  In 
this  case,  the  optimization  problem  whose  solution  speci¬ 
fies  the  sound  speed  estimate  can  be  re-expressed  in  a 
stochastic  framework  as 

c  =  argrninEw [1  ||  gw  -  MHcsw  ||2}  +  /311(c),  (15) 

which  we  refer  to  as  the  WISE  method.  An  implemen¬ 
tation  of  the  WISE  method  that  utilizes  the  stochastic 
gradient  descent  algorithm  is  summarized  in  Algorithm  2. 

In  Algorithm  2,  the  wave  equation  solver  needs  to  be 
run  one  time  in  each  of  Lines  5  and  6.  In  the  line  search 
to  determine  the  step  size  in  Line  8,  the  wave  equation 
solver  needs  to  be  run  at  least  once,  but  in  general  will 
require  a  small  number  of  additional  runs,  just  as  in  Al¬ 
gorithm  1.  Accordingly,  the  lower  bound  on  the  number 
of  required  wave  equation  solver  runs  per  iteration  is  3,  as 
opposed  to  (3 M)  for  the  conventional  sequential  waveform 
inversion  method  described  by  Algorithm  1.  As  demon¬ 
strated  in  geophysics  applications  [29],  [31],  [41]  and  the 
breast  imaging  studies  below,  the  WISE  method  provides 
a  substantial  reduction  in  reconstruction  times  over  use 
of  the  standard  sequential  waveform  inversion  method.  In 
Line  7  of  Algorithm  2,  Jw  can  be  calculated  analogously 
to  (10)  as 

rrwi  _  1  Wfnwl  [VW]nL+l-l  ~  %[VW]nL+l  +  [VW]nL+l+l 

LJ  L~rj3  2Jq  J  nL+(L-l)  Tt  ’ 

lc\n  i= i  ^ 

(16) 

where  pw  =  Hcsw  and  qw  =  Hcrw  with  rw  e  RNL  calcu¬ 
lated  by 


A.  Measurement  Geometry 

A  circular  measurement  geometry  was  chosen  to  emu¬ 
late  a  previously  reported  USCT  breast  imaging  system 
[10],  [23],  [43].  As  depicted  in  Fig.  1,  256  ultrasonic  trans¬ 
ducers  were  uniformly  distributed  on  a  ring  of  radius 
110  mm.  The  generation  of  one  USCT  data  set  consisted 
of  M  =  256  sequential  data  acquisitions.  In  each  data  ac¬ 
quisition,  one  emitter  produced  an  acoustic  pulse.  The 
acoustic  pulse  was  numerically  propagated  through  the 
breast  phantom  and  the  resulting  wavefield  data  were  re¬ 
corded  by  all  transducers  in  the  array  as  described  below. 
Note  that  the  location  of  the  emitter  in  every  data  ac¬ 
quisition  was  different  from  those  in  other  acquisitions, 
whereas  the  locations  of  receivers  were  identical  for  all 
acquisitions. 

B.  Numerical  Breast  Phantom 

A  numerical  breast  phantom  of  diameter  98  mm  was 
employed.  The  phantom  was  composed  of  8  structures 
representing  adipose  tissues,  parenchymal  breast  tissues, 
cysts,  benign  tumors,  and  malignant  tumors,  as  shown  in 
Fig.  2.  For  simplicity,  the  acoustic  attenuation  of  all  tis¬ 
sues  was  described  by  a  power  law  with  a  fixed  exponent 
y  =  1.5  [44].  The  corresponding  sound  speed  and  the  at¬ 
tenuation  slope  values  are  listed  in  Table  I  [44]— [46] .  Both 
the  sound  speed  and  the  attenuation  slope  distributions 
in  Fig.  2  were  sampled  on  a  uniform  Cartesian  grid  with 
spacing  As  =  0.25  mm.  The  finest  structure  [indexed  by  7 
in  Fig.  2(a)]  was  of  diameter  3.75  mm. 

C.  Simulation  of  the  Measurement  Data 

1 )  First-Order  Numerical  Wave  Equation  Solver:  Acous¬ 
tic  wave  propagation  in  acoustically  absorbing  media  was 
modeled  by  three  coupled  first-order  partial  differential 
equations  [47]: 
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(a) 


Fig.  2.  (a)  Sound  speed  map  (mm-ps  x)  and  (b)  acoustic  attenuation 
slope  map  [dB-(MHz)~Vcm_1]  Qf  the  numerical  breast  phantom. 


TABLE  I.  Parameters  of  the  Numerical  Breast 
Phantom  [44]— [46] . 


Structure 

index 

Tissue  type 

Sound  speed 
(mm-ps-1) 

Slope  of  attenuation 
[dB-(MHz)-vCm-1] 

0 

Adipose 

1.47 

0.60 

1 

Parenchyma 

1.51 

0.75 

2 

Benign  tumor 

1.47 

0.60 

3 

Benign  tumor 

1.47 

0.60 

4 

Cyst 

1.53 

0.00217 

5 

Malignant  tumor 

1.565 

0.57 

6 

Malignant  tumor 

1.565 

0.57 

7 

Malignant  tumor 

1.57 

0.57 

^  u(r,i)  =  -Vp{r,t),  (18a) 

if p{r,t )  =  -V  •  u(r,i)  +  4irfdt's(r,t'),  (18b) 

p(r,t)  = 

c2(r)[l  +  r(r)A(_v2)^-i  +  v( r)(-V^y+^]p{r,t), 

(18c) 

where  u(r,t),  p(r,£),  and  p{ r)  denote  the  acoustic  particle 
velocity,  the  acoustic  pressure,  and  the  acoustic  density, 
respectively.  The  functions  r(r)  and  77 (r)  describe  acoustic 
absorption  and  dispersion  during  the  wave  propagation 
[47]: 

r(r)  =  -2a0(r)c0(r)!'"1,  r?(r)  =  2a0(r)c0(r)!' tan(7rj//2), 

(19) 

where  a'o(r)  and  y  are  the  attenuation  slope  and  the  power 
law  exponent,  respectively.  When  the  medium  is  assumed 
to  be  lossless  [i.e.,  a'o(r)  =  0],  it  can  be  shown  that  (18)  is 
equivalent  to  (1). 

Based  on  (18),  a  pseudospectral  k-space  method  was 
employed  to  simulate  acoustic  pressure  data  [36],  [47]. 
This  method  was  implemented  by  use  of  a  first-order  nu¬ 
merical  scheme  on  GPU  hardware.  The  calculation  do¬ 
main  was  of  size  512  x  512  mm2,  sampled  on  a  2048  x 
2048  uniform  Cartesian  grid  of  spacing  As  =  0.25  mm.  A 
nearest-neighbor  interpolation  was  employed  to  place  all 


Fig.  3.  (a)  Normalized  temporal  profile  and  (b)  amplitude  spectrum  of 
the  excitation  pulse  employed  in  the  computer  simulation  studies.  The 
dashed  line  in  (b)  marks  the  center  frequency  of  excitation  pulse  at 
0.82  MHz. 


transducers  on  the  grid  points.  On  a  platform  consisting 
of  dual  quad-core  CPUs  with  a  3.30  GHz  clock  speed, 
64  GB  of  random-access  memory,  and  a  single  NVIDIA 
Tesla  K20  GPU,  the  first-order  pseudospectral  k-space 
method  required  approximately  108  s  to  complete  one  for¬ 
ward  simulation. 


2)  Acoustic  Excitation  Pulse:  The  excitation  pulse  em¬ 
ployed  in  this  study  was  assumed  to  be  spatially  localized 
at  the  emitter  location  while  temporally  it  was  a  fc  = 
0.8  MHz  sinusoidal  function  tapered  by  a  Gaussian  kernel 
with  standard  deviation  a  =  0.5  ps,  i.e., 


*».(*■•  0  = 


(t-tc)2} 

exp 

.0, 

2<t2  J 

sin(27r  fct), 


at  the  rath  emitter  location, 
otherwise, 

(20) 


where  tc  =  3.2  ps.  The  temporal  profile  and  the  amplitude 
frequency  spectrum  of  the  excitation  pulse  are  plotted  in 
Fig.  3(a)  and  (b),  respectively.  The  excitation  pulse  con¬ 
tained  approximately  3  cycles. 


3)  Generation  of  Non- Attenuated  and  Attenuated  Noise- 
Free  Data:  For  every  data  acquisition  (indexed  by  ra),  the 
first-order  pseudospectral  k-space  method  was  run  for 
3600  time  steps  with  a  time  interval  At  =  0.05  ps  (corre¬ 
sponding  to  a  20  MHz  sampling  rate).  Downsampling  the 
recorded  data  by  taking  every  other  time  sample  resulted 
in  a  data  vector  gm  that  was  effectively  sampled  at 
10  MHz  and  was  of  dimensions  ML  with  M  =  256  and  L 
=  1800.  The  data  vector  at  the  zeroth  data  acquisition,  g0, 
is  displayed  as  a  2-D  image  in  Fig.  4(a).  This  undersam¬ 
pling  procedure  was  introduced  to  avoid  inverse  crime  [48] 
so  that  the  data  generation  and  the  image  reconstruction 
employed  different  numerical  discretization  schemes.  Re¬ 
peating  the  calculation  for  ra  =  0,  1,  ...,  255,  we  obtained 
a  collection  {gm}  of  data  vectors  that  together  represented 
one  complete  data  set.  Using  the  absorption  phantom  de¬ 
scribed  in  Section  IV-B,  a  complete  attenuated  data  set 
was  computed.  An  idealized,  non- attenuated  data  set  was 
also  computed  by  setting  a'o(r)  =  0. 
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(a)  (b)  (c) 

Fig.  4.  Computer-simulated  (a)  noise-free  and  (b)  noisy  data  vectors  at  the  zeroth  data  acquisition,  (c)  Profiles  of  the  pressure  received  by  the  128th 
transducer.  The  grayscale  window  for  (a)  and  (b)  is  [—45,  0]  dB. 


4 )  Generation  of  Incomplete  Data:  An  incomplete  data 
set  in  this  study  corresponds  to  one  in  which  only  Nrec 
receivers  located  on  the  opposite  side  of  the  emitter  re¬ 
cord  the  pressure  wavefield,  with  NTec  <  M.  Taking  the 
zeroth  data  acquisition  as  an  example  (Fig.  1),  only  Nrec 
=  100  receivers,  indexed  from  78  to  177,  record  the  wave- 
field,  whereas  other  receivers  record  either  unreliable  or 
no  measurements.  Incomplete  data  sets  formed  in  this  way 
can  emulate  two  practical  scenarios:  (1)  Signals  recorded 
by  receivers  near  the  emitter  are  unreliable  and  therefore 
discarded  [23],  and  (2)  an  arc-shaped  transducer  array  is 
employed  that  rotates  with  the  emitter  [13],  [14],  [49]. 

Specifically,  incomplete  data  sets  were  generated  as 

.  .  m  =  0,1,  •••,  M  —  1 

[g mCP  ]„■*£  +  ;  =  [g mbjOi  +  i’  for  =  0  1  •••  VreC  -  l’ 

(21) 

where  g™cpl  is  the  incomplete  rath  data  acquisition,  which 
is  of  dimensions  ArecL,  with  Nrec  <  M.  The  index  map 
Jm  :  {0,1,  *  •  • , Nrec  1 }  i  ^  M^od  is  defined  as 

I  M  -  Nrec\ 

Jm{nrec )  =  [ra  +  nrec  + - g - j  mod  M ,  (22) 

where  (ra'  mod  M)  calculates  the  remainder  of  m'  divided 
by  M,  and  the  index  set  M^od  collects  indices  of  transduc¬ 
ers  that  reliably  record  data  at  the  rath  data  acquisition 
and  is  defined  as 


{k  mod  M  \k  e[m  +  (M  -  Arec)/2,ra  +  (M  +  Arec)/2)}. 

(23) 

Here,  for  simplicity,  we  assume  that  M  and  N*ec  are  even 
numbers.  In  this  study,  we  empirically  set  A4ec  =  100  so 
that  the  object  can  be  fully  covered  by  the  fan  region  as 
shown  in  Fig.  1. 

5)  Generation  of  Noisy  Data:  An  additive  Gaussian 
white  noise  model  was  employed  to  simulate  electronic 
measurement  noise  as 


|m  =  gm  +  n,  (24) 

where  gm  and  n  are  the  noisy  data  vector  and  the  Gauss¬ 
ian  white  noise  vector,  respectively.  In  this  study,  the 
maximum  value  of  the  pressure  received  by  the  128th 
transducer  at  the  zeroth  data  acquisition  with  a  homoge¬ 
neous  medium  (water  tank)  was  chosen  as  a  reference  sig¬ 
nal  amplitude.  The  noise  standard  deviation  was  set  to  be 
5%  of  this  value.  An  example  of  a  simulated  noiseless  and 
noisy  data  acquisition  is  shown  Fig.  4. 

D.  Image  Reconstruction 

1 )  Second- Order  Pseudo  spectral  k-  Space  Method:  In  the 
reconstruction  methods  described  below,  the  action  of  the 
operator  Hc  (5)  was  computed  by  solving  (1)  by  use  of 
a  second-order  psendospectral  k-space  method.  This  was 
implemented  using  GPUs.  The  calculation  domain  was  of 
size  512  x  512  mm2,  sampled  on  a  1024  x  1024  uniform 
Cartesian  grid  of  spacing  As  =  0.5  mm  for  reconstruction. 
On  a  platform  consisting  of  dual  octa-core  CPUs  with  a 
2.00  GHz  clock  speed,  125  GB  RAM,  and  a  single  NVID¬ 
IA  Tesla  K20C  GPU,  the  second-order  k-space  method 
required  approximately  7  s  to  complete  one  forward  simu¬ 
lation. 

2)  Sequential  Waveform  Inversion:  To  serve  as  a  refer¬ 
ence  for  the  WISE  method,  we  implemented  the  sequen¬ 
tial  waveform  inversion  method  described  in  Algorithm 
1.  A  uniform  sound  speed  distribution  was  employed  as 
the  initial  guess,  which  corresponded  to  the  known  back¬ 
ground  value  of  1.5  mm/ps.  The  object  was  contained 
in  a  square  region  of  interest  (ROI)  of  dimension  128  x 
128  mm2  (Fig.  1),  which  was  covered  by  256  x  256  pixels. 

3)  WISE  Method:  We  implemented  the  WISE  method 
by  use  of  Algorithm  2.  Two  types  of  penalties  were  em¬ 
ployed  in  this  study:  a  quadratic  penalty  expressed  as 

^Q(c)  =  £B«CW<  -  [cWh-i)8 

3  i  (25) 

+  (MjiV.+i  _  [<%— l)iVI+i)2}, 
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Vec)  = 


M  —  Niec  ,rM-Niec  rpr  M  +  Niec 

- 2 - ’  lf - 2 - -m  ~m< - 2 - 

rec  M  +  iVrec  -M  -  Nlec  ^  rec  —M  +  NTec 

m  -  m  H - = - ,  if - x - <  m  -  m  < - s - . 


(28) 


where  Nx  and  Ny  denote  the  number  of  grid  points  along 
the  x  and  y  directions,  respectively,  and  a  total  variation 
(TV)  penalty,  defined  as  [52],  [53] 

nTV(c)  = 

XTa/£  +  (WiJVj  +  i  _  +  +  (WjA^  +  i  -  My-ijjv^+i)2 , 

3  i 

(26) 


strategy  is  based  on  the  assumption  that  the  back-scatter 
from  breast  tissue  in  an  appropriately  sound  speed- 
matched  water  bath  is  weak.  This  assumption  suggests 
that  the  missing  measurements  can  be  replaced  by  the 
corresponding  pressure  data  that  would  have  been  pro¬ 
duced  in  the  absence  of  the  object. 

The  second,  more  crude,  data  completion  strategy  was 
to  simply  fill  the  missing  data  with  zeros,  i.e., 


where  e  is  a  small  number  introduced  to  avoid  dividing 
by  0  in  the  gradient  calculation.  In  this  study,  we  empiri¬ 
cally  selected  e  =  10~12.  This  value  was  fixed  because  we 
observed  that  it  had  a  minor  impact  on  the  reconstructed 
images  compared  with  the  impact  of  j3.  The  use  of  this  pa¬ 
rameter  can  be  avoided  when  advanced  optimization  algo¬ 
rithms  are  employed  [54],  [55].  As  in  the  sequential  wave¬ 
form  inversion  case,  it  was  assumed  that  the  background 
sound  speed  was  known  and  the  object  was  contained  in  a 
square  ROI  of  dimension  128  x  128  mm2  (Fig.  1),  which 
corresponded  to  256  x  256  pixels.  The  regularization  pa¬ 
rameters  corresponding  to  the  quadratic  penalty  and  the 
TV  penalty  will  be  denoted  by  f3®  and  /3TV,  respectively. 
Optimal  regularization  parameter  values  should  ultimate¬ 
ly  be  identified  by  use  of  task-based  measures  of  image 
quality  [37].  In  this  preliminary  study,  we  investigated  the 
impact  of  (3®  and  /3TV  on  the  reconstructed  images  by 
sweeping  their  values  over  a  wide  range. 

4)  Reconstruction  From  Incomplete  Data:  Because  the 
WISE  method  requires  Mm  to  be  identical  for  all  ra’s, 
image  reconstruction  from  incomplete  data  remains  chal¬ 
lenging  [30],  [33],  [42].  In  this  study,  two  data  completion 
strategies  were  investigated  [30],  [33],  [42]  to  synthesize  a 
complete  data  set,  to  which  the  WISE  method  could  be 
effectively  applied. 

One  strategy  was  to  fill  the  missing  data  with  pressure 
corresponding  to  a  homogeneous  medium  as 

r  combH,  _  ( \^\,WL+V  if  *  if 

LSm  J m™L  +  l  ~  1  r  h ,  , , 

[[gjm,eci+p  otherwise, 

(27) 

for  mTec  =  0,  1, M -  1,  where  g*  G  RML,g^cpl  €  R1*™1, 
and  g“mbH  G  RML  denote  the  computer-simulated  (with  a 
homogeneous  medium),  the  measured  incomplete,  and  the 
combined  complete  data  vectors  at  the  rath  data  acquisi¬ 
tion,  respectively.  The  mapping 

Jml  :  M^od  i  ^  {0, 1 ,  -  -  - ,  Arec  —  1}  denotes  the  inverse  op¬ 
erator  of  Jm  as  in  (28),  see  above.  This  data  completion 


r  combO -I 

18m  J  miecL+l 


^mapli 

o, 


if  mrec  G  MS°od 
otherwise, 

(29) 


where  g^mb°  denotes  the  data  completed  with  the  second 
strategy. 


5)  Bent- Ray  Image  Reconstruction:  A  bent-ray  method 
was  also  employed  to  reconstruct  images.  Details  regard¬ 
ing  the  TOF  estimation  and  algorithm  implementation 
are  provided  in  Appendix  C. 


V.  Computer  Simulation  Results 
A.  Images  Reconstructed  From  Idealized  Data 

The  images  reconstructed  from  the  noise-free,  non-at- 
tennated  data  by  use  of  the  WISE  method  with  199  itera¬ 
tions  and  the  sequential  waveform  inversion  method  with 
43  iterations  are  shown  in  Fig.  5(a)  and  (b).  As  expected 
[23] ,  [56] ,  both  images  are  more  accurate  and  possess  higher 
spatial  resolution  than  the  one  reconstructed  by  use  of  the 
bent-ray  reconstruction  algorithm  displayed  in  Fig.  5(c). 
Profiles  through  the  reconstructed  images  are  displayed 
in  Fig.  6.  The  images  shown  in  Fig.  5(a)  and  (b)  possess 
similar  accuracies  as  measured  by  their  root-mean-sqnare 
errors  (RMSEs),  namely,  1.08  x  10-3  for  the  former  and 
1.19  x  10-3  for  the  latter.  The  RMSE  was  computed  as 
the  Euclidean  distance  between  the  reconstructed  image 
and  the  sound  speed  phantom  vector  c,  averaged  by  the 
256  x  256  pixels  of  the  ROI  sketched  in  Fig.  1.  However, 
the  reconstruction  of  Fig.  5(a)  required  only  about  1.7% 
of  the  computational  time  required  to  reconstruct  Fig. 
5(b),  namely,  1.4  hours  for  the  former  and  81.4  hours 
for  the  latter.  This  is  because  the  WISE  method  required 
only  1018  wave  equation  solver  runs,  which  is  significantly 
less  than  the  57088  wave  equation  solver  runs  required  by 
the  sequential  waveform  inversion  method.  With  a  similar 
number  of  wave  equation  solver  runs  (e.g.,  1024),  one  can 
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(C)  (d) 

Fig.  5.  Images  reconstructed  by  use  of  (a)  the  WISE  method  after  the 
199th  iteration  (1018  runs  of  the  wave  equation  solver),  (b)  the  sequen¬ 
tial  waveform  inversion  algorithm  after  the  43rd  iteration  (57088  runs 
of  the  wave  equation  solver),  (c)  the  bent-ray  model-based  sound  speed 
reconstruction  method,  and  (d)  the  sequential  waveform  inversion  algo¬ 
rithm  after  the  first  iteration  (1024  runs  of  the  wave  equation  solver) 
from  the  noise- free  non-attenuated  data.  The  grayscale  window  is  [1.46, 
1.58]  mm/fxs. 


complete  only  a  single  algorithm  iteration  by  use  of  the 
sequential  waveform  inversion  method.  The  corresponding 
image,  shown  in  Fig.  5(d),  lacks  quantitative  accuracy  as 
well  as  qualitative  value  for  identifying  features.  The  re¬ 
sults  suggest  that  the  WISE  method  maintains  the  advan¬ 
tages  of  the  sequential  waveform  inversion  method  while 
significantly  reducing  the  computational  time. 

B.  Convergence  of  the  WISE  Method 

Images  reconstructed  from  noise-free,  non-attenuated, 
data  by  use  of  the  WISE  method  contain  radial  streak 
artifacts  when  the  algorithm  iteration  number  is  less  than 
100,  as  shown  in  Figs.  7(a)  to  (c).  Profiles  through  these 
images  are  displayed  in  Fig.  8.  The  streak  artifacts  are 
likely  caused  by  crosstalk  introduced  during  the  source 
encoding  procedure  [31],  [41].  However,  these  artifacts  are 
effectively  mitigated  after  more  iterations  as  demonstrat¬ 
ed  by  the  image  reconstructed  after  the  199th  iteration 
in  Fig.  5(a)  and  its  profile  in  Fig.  6.  The  quantitative  ac¬ 
curacy  of  the  reconstructed  images  is  improved  with  more 
iterations  as  shown  in  Fig.  8. 

Fig.  9(a)  reveals  that  the  WISE  method  requires  a  larg¬ 
er  number  of  algorithm  iterations  than  does  the  sequential 
waveform  inversion  method  to  achieve  the  same  RMSE. 


Fig.  6.  Profiles  at  y  —  6.5  mm  of  the  images  reconstructed  by  use  of  the 
bent-ray  TOF  image  reconstruction  method  and  the  WISE  method  from 
the  noise-free  non-attenuated  data. 


The  RMSE  of  the  images  reconstructed  by  use  of  the 
WISE  method  appears  to  oscillate  around  1.0  x  10-3  after 
the  first  100  iterations,  whereas  the  sequential  waveform 
inversion  method  can  achieve  a  lower  RMSE.  However,  as 
shown  previously  in  Fig.  5(a)  and  the  corresponding  pro¬ 
file  in  Fig.  6,  after  additional  iterations  the  image  recon¬ 
structed  by  use  of  the  WISE  method  achieves  a  high  ac¬ 
curacy.  Moreover,  to  achieve  the  same  accuracy  as  the 
sequential  waveform  inversion  method,  the  WISE  method 
requires  a  computation  time  that  is  reduced  by  approxi¬ 
mately  2  orders  of  magnitude,  as  suggested  by  Fig.  9(b). 
We  also  plotted  the  cost  function  value  against  the  num¬ 
ber  of  iterations  in  Fig.  9(c).  Note  that  for  the  WISE 
method,  the  cost  function  value  was  approximated  by  the 

current  realization  of  ||g  —  MHcsw||  /2.  These  plots  sug¬ 

gest  that,  in  this  particular  case,  the  WISE  method  ap¬ 
pears  to  approximately  converge  after  200  iterations.  For 
example,  the  images  reconstructed  after  199  [Fig.  5(a)] 
and  250  [Fig.  7(d)]  iterations  are  nearly  identical. 

C.  Images  Reconstructed  From  Non- Attenuated 
Data  Containing  Noise 

Images  reconstructed  by  use  of  the  WISE  method  with 
a  quadratic  penalty  and  the  WISE  method  with  a  TV 
penalty  from  noisy,  non-attenuated  data  are  presented  in 
Fig.  10.  All  images  were  obtained  after  1024  algorithm 
iterations.  The  WISE  method  with  a  quadratic  penalty  ef¬ 
fectively  mitigates  image  noise  as  shown  in  Figs.  10(a)  to 
(c),  at  the  expense  of  image  resolution,  as  expected.  Fig. 
10(d)  shows  an  image  reconstructed  by  use  of  the  WISE 
method  with  a  TV  penalty.  The  image  appears  to  possess 
a  similar  resolution  but  a  lower  noise  level  than  the  image 
in  Fig.  10(b)  that  was  reconstructed  by  use  of  the  WISE 
method  with  a  quadratic  penalty.  We  also  compared  the 
convergence  rates  of  the  WISE  method  and  the  sequential 
waveform  inversion  methods  when  both  utilize  a  TV  pen¬ 
alty  and  the  same  regularization  parameter.  As  shown  in 
Fig.  11,  the  convergence  properties  of  the  penalized  meth- 
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(C)  (d) 

Fig.  7.  Images  reconstructed  by  use  of  the  WISE  method  after  (a)  the 
20th,  (b)  the  50th,  (c)  the  100th,  and  (d)  the  250th  iteration  from  the 
noise- free,  non- attenuated  data  set.  The  grayscale  window  is  [1.46, 
1.58]  mm/qs. 


ods  follow  similar  trends  as  the  un-penalized  methods, 
which  were  discussed  above  and  shown  in  Fig.  9.  Even 
though  it  required  a  larger  number  of  algorithm  iterations, 
the  WISE  method  reduced  the  computation  time  by  ap¬ 
proximately  2  orders  of  magnitude  as  compared  with  the 
sequential  waveform  inversion  method. 

D.  Images  Reconstructed  From  Acoustically 
Attenuated  Data 

Our  current  implementation  of  the  WISE  method  as¬ 
sumes  an  absorption-free  acoustic  medium.  This  assump¬ 


x  [mm] 

Fig.  8.  Profiles  of  the  images  reconstructed  by  use  of  the  WISE  method 
from  the  noise-free  non- attenuated  data  after  different  numbers  of  itera¬ 
tions. 


tion  can  be  strongly  violated  in  practice.  To  investigate 
the  robustness  of  the  WISE  method  to  model  errors  as¬ 
sociated  with  ignoring  medium  acoustic  absorption,  we 
applied  the  algorithm  to  the  acoustically  attenuated  data 
that  were  produced  as  described  in  Section  IV-C.  As 
shown  in  Fig.  12,  when  acoustic  absorption  is  considered, 
the  amplitude  of  the  measured  pressure  is  attenuated  by 
approximately  a  factor  of  2.  The  wavefront  [Fig.  12(a)] 
remains  very  similar  to  that  when  medium  absorption  is 
ignored  [Fig.  4(a)].  Medium  absorption  has  the  largest  im¬ 
pact  on  the  pressure  data  received  by  transducers  located 
opposite  the  emitter  as  shown  in  Fig.  12(b).  The  shape  of 
the  pulse  profile  remains  very  similar  as  shown  in  Figs. 
12(c)  and  (d),  suggesting  that  waveform  dispersion  may 
be  less  critical  than  amplitude  attenuation  in  image  recon¬ 
struction  for  this  phantom. 

Images  reconstructed  by  use  of  the  WISE  method  with 
a  TV  penalty  from  noise-free  and  noisy  attenuated  data 
are  shown  in  Figs.  13(a)  and  (b).  Image  profiles  are  shown 
in  Fig.  13(c).  Although  these  images  contain  certain  arti¬ 
facts  that  were  not  produced  in  the  idealized  data  studies, 
most  object  structures  remain  readily  identified.  These  re¬ 
sults  suggest  that  the  WISE  method  with  a  TV  penalty 
can  tolerate  data  inconsistencies  associated  with  neglect- 


Fig.  9.  Plots  of  the  root  mean  square  errors  (RMSEs)  of  the  images  reconstructed  from  the  noise-free  data  versus  (a)  the  number  of  iterations  and 
(b)  the  number  of  wave  equation  solver  runs,  (c)  Plots  of  the  cost  function  value  versus  the  number  of  iterations. 
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(c)  (d) 

Fig.  10.  Images  reconstructed  from  non- attenuated  data  contaminated 
with  Gaussian  random  noise.  Images  (a)  to  (c)  were  reconstructed  by  use 
of  the  WISE  method  with  a  quadratic  penalty  with  {D  =  1.0  x  10-3, 
1.0  x  10-2,  and  1.0  x  10_1,  respectively.  Image  (d)  was  reconstructed  by 
use  of  the  WISE  method  with  a  TV  penalty  with  /dTV  =  5.0  x  10-4.  The 
insert  in  the  up  right  corner  of  each  image  is  the  zoomed-in  image  of  the 
dashed  black  box,  which  contains  35  x  35  pixels  (17.5  x  17.5  mm2).  The 
grayscale  window  is  [1.46,  1.58]  mm/ps. 


ing  acoustic  attenuation  in  the  imaging  model,  at  least  to 
a  certain  level  with  regard  to  feature  detection  tasks. 

E.  Images  Reconstructed  From  Idealized  Incomplete  Data 

The  wavefront  of  the  noise-  and  attenuation- free  pres¬ 
sure  wavefield  when  the  object  is  absent  [Fig.  14(a)]  ap¬ 


pears  to  be  very  similar  to  that  when  the  object  is  present 
[Fig.  4(a)].  As  expected,  the  largest  differences  are  seen  in 
the  signals  received  by  the  transducers  located  opposite 
of  the  emitter,  as  shown  in  Fig.  14(b).  As  seen  in  Fig. 
14(c),  the  time  traces  received  by  the  40th  transducer  are 
nearly  identical  when  object  is  present  and  absent.  This 
is  because  the  back-scattered  wavefield  is  weak  for  breast 
imaging  applications.  These  results  establish  the  potential 
efficacy  of  the  data  completion  strategy  of  filling  the  miss¬ 
ing  data  with  the  pressure  data  corresponding  to  a  water 
bath. 

The  image  reconstructed  from  the  measurements  com¬ 
pleted  with  pressure  data  corresponding  to  a  water  bath 
is  shown  in  Fig.  15(a).  As  revealed  by  the  profile  in  Fig. 
15(c),  this  image  is  highly  accurate.  Alternatively,  the 
image  reconstructed  from  the  data  completed  with  zeros 
contains  strong  artifacts  as  shown  in  Fig.  15(b).  These 
results  suggest  that  the  WISE  method  can  be  adapted  to 
reconstruct  images  from  incomplete  data,  which  is  par¬ 
ticularly  useful  for  emerging  laser-induced  USCT  imaging 
systems  [13]— [15] . 


VI.  Experimental  Validation 
A.  Data  Acquisition 

Experimental  data  recorded  by  use  of  the  SoftVue 
USCT  scanner  [57]  were  utilized  to  further  validate  the 
WISE  method.  The  scanner  contained  a  ring-shaped  array 
of  radius  110  mm  that  was  populated  with  2048  trans¬ 
ducer  elements.  Each  element  had  a  center  frequency  of 
2.75  MHz,  a  pitch  of  0.34  mm,  and  was  elevationally  fo¬ 
cused  to  isolate  a  3-mm-thick  slice  of  the  to-be-imaged 
object.  The  transducer  array  was  mounted  in  a  water  tank 
and  could  be  translated  with  a  motorized  gantry  in  the 
vertical  direction.  See  [57]  for  additional  details  regarding 
the  system. 

The  breast  phantom  was  built  by  Dr.  Ernie  Madsen 
from  the  University  of  Wisconsin  and  provides  tissue- 


Fig.  11.  Plots  of  the  root- mean-square  errors  (RMSEs)  of  the  images  reconstructed  from  the  noisy  data  versus  (a)  the  number  of  iterations  and  (b) 
the  number  of  wave  equation  solver  runs,  (c)  Plots  of  the  cost  function  value  versus  the  number  of  iterations.  Both  the  WISE  and  the  sequential 
waveform  inversion  methods  employed  a  TV  penalty  with  /7rv  =  5.0  x  10-4. 
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Fig.  12.  (a)  Computer-simulated  noise- free  attenuated  pressure  of  the 
zeroth  data  acquisition,  (b)  The  difference  between  the  attenuated  pres¬ 
sure  data  and  the  non- attenuated  pressure  data,  (c)  The  temporal  pro¬ 
files  and  (d)  the  amplitude  spectra  of  the  pressure  received  by  the  128th 
transducer.  The  grayscale  window  for  (a)  and  (b)  is  [—45,  0]  dB. 


equivalent  characteristics  of  highly  scattering,  predomi¬ 
nantly  parenchymal  breast  tissue.  The  phantom  mimics 
the  presence  of  benign  and  cancerous  masses  embedded  in 
glandular  tissue,  including  a  subcutaneous  fat  layer.  Fig. 
16  displays  a  schematic  of  one  slice  through  the  phantom. 
The  diameter  of  the  inclusions  is  approximately  12  mm. 
Table  II  presents  the  known  acoustic  properties  of  the 
phantom. 

During  data  acquisition,  the  breast  phantom  was 
placed  near  the  center  of  the  ring-shaped  transducer  ar¬ 
ray  so  that  the  distance  between  the  phantom  and  each 
transducer  was  approximately  the  same.  While  scanning 
each  slice,  every  other  transducer  element  sequentially 
emits  fan  beam  ultrasound  signals  toward  the  opposite 
side  of  the  ring.  The  forward-scattered  and  back-scattered 
ultrasound  signals  are  subsequently  recorded  by  the  same 
transducer  elements.  The  received  waveform  was  sampled 
at  a  rate  of  12  MHz.  The  1024  data  acquisitions  required 
approximately  20  s  in  total.  A  calibration  data  set  was 
also  acquired  in  which  the  phantom  object  was  absent. 


(a) 


(b) 


(c) 


Fig.  13.  (a)  Image  reconstructed  by  use  of  the  WISE  method  from  the  noise- free  attenuated  data,  (b)  Image  reconstructed  by  use  of  the  WISE 
method  with  a  TV  penalty  with  dTV  =  5.0  x  10-4,  from  the  noisy  attenuated  data.  The  grayscale  window  is  [1.46,  1.58]  mm/ps.  (c)  Profiles  at  y 
=  6.5  mm  of  the  reconstructed  images. 


Fig.  14.  (a)  Computer-simulated  noise-free  non-attenuated  pressure  data  when  the  object  is  absent,  (b)  The  difference  between  the  pressure  data 
when  object  is  present  and  the  pressure  data  when  the  object  is  absent,  (c)  Profiles  of  the  pressure  received  by  the  40th  transducer.  The  grayscale 
window  for  (a)  and  (b)  is  [—45,  0]  dB. 
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(a)  (b)  (c) 

Fig.  15.  Images  reconstructed  by  use  of  the  WISE  method  from  noise-free  combined  data  that  are  completed  (a)  with  computer-simulated  pressure 
corresponding  to  a  homogeneous  medium  and  (b)  with  zeros.  The  grayscale  window  is  [1.46,  1.58]  mm/ps.  (c)  Profiles  at  y  =  6.5  mm  of  the  images 
reconstructed  by  use  of  the  WISE  method  from  the  two  combined  data  sets. 


B.  Data  Pre-Processing 

Forty-eight  bad  channels  were  manually  identified  by 
visual  inspection.  After  discarding  these,  the  data  set  con¬ 
tained  M  =  976  acquisitions.  Each  acquisition  contained 
A rec  =  976  time  traces.  Each  time  trace  contained  L  = 
2112  time  samples.  The  976  good  channels  were  indexed 
from  0  to  975.  The  corresponding  data  acquisitions  were 
indexed  in  the  same  way.  A  Hann-window  low-pass  filter 
with  a  cutoff  frequency  of  4  MHz  was  applied  to  every 
time  trace  in  both  the  calibration  and  the  measurement 
data.  This  data  filtering  was  implemented  to  mitigate  nu¬ 
merical  errors  that  could  be  introduced  by  our  second- 
order  wave  equation  solver. 

C.  Estimation  of  Excitation  Pulse 

The  shape  of  the  excitation  pulse  was  estimated  as  the 
time  trace  of  the  calibration  data  (after  pre-processing) 
received  by  the  488th  receiver  at  the  zeroth  data  acquisi¬ 
tion.  Note  that  the  488th  receiver  was  approximated  lo¬ 
cated  on  the  axis  of  the  zeroth  emitter,  thus  the  received 
pulse  was  minimally  affected  by  the  finite  aperture  size 


Fibroadenoma 


Fig.  16.  Schematic  of  the  breast  phantom  employed  in  the  experimental 
study. 


effect  of  the  transducers.  Because  our  calibration  data  and 
measurement  data  were  acquired  using  different  electronic 
amplifier  gains,  the  amplitude  of  the  excitation  pulse  was 
estimated  from  the  measurement  data.  More  specifically, 
we  simulated  the  zeroth  data  acquisition  using  the  sec¬ 
ond-order  pseudospectral  k-space  method  and  compared 
the  simulated  time  trace  received  by  the  300th  receiver 
with  the  corresponding  measured  time  trace  (after  pre¬ 
processing).  The  ratio  between  the  maximum  values  of 
these  two  traces  was  used  to  scale  the  excitation  pulse 
shape.  We  selected  the  300th  receiver  because  it  resided 
out  of  the  fan  region  indicated  in  Fig.  1;  its  received  sig¬ 
nals  were  unlikely  to  be  strongly  affected  by  the  presence 
of  the  object.  The  estimated  excitation  pulse  and  its  am¬ 
plitude  spectrum  are  displayed  in  Fig.  17.  Note  that  the 
experimental  excitation  pulse  contained  higher  frequency 
components  than  did  the  computer-simulated  excitation 
pulse  shown  in  Fig.  3. 

D.  Synthesis  of  Combined  Data 

As  discussed  in  Section  IV-C-4,  signals  received  by  re¬ 
ceivers  located  near  the  emitter  can  be  unreliable  [23]. 
Our  experimental  data,  as  shown  in  Fig.  18(a),  contained 
noise-like  measurements  for  the  receivers  indexed  from  0 
to  200,  and  from  955  to  975,  in  the  case  where  the  zeroth 
transducer  functioned  as  the  emitter.  Also,  our  point-like 
transducer  assumption  introduces  larger  model  mismatch- 


TABLE  II.  Parameters  of  the  Experimental  Breast  Phantom. 


Material 

Sound  speed 
(mm-ps 1) 

Attenuation 
coefficient  at  2.5 

MHz  (dB/cm) 

Fat 

1.467 

0.48 

Parenchymal  tissue 

1.552 

0.89 

Cancer 

1.563 

1.20 

Fibroadenoma 

1.552 

0.52 

Gelatin  cyst 

1.585 

0.16 
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(a)  (b) 

Fig.  17.  (a)  Normalized  temporal  profile  and  (b)  amplitude  spectrum  of 
the  excitation  pulse  employed  in  the  experimental  studies.  The  dashed 
line  in  (b)  marks  the  center  frequency  of  excitation  pulse  at  2.09  MHz. 

es  for  the  receivers  located  near  the  emitter.  As  shown 
in  Figs.  18(c)  and  (d),  even  though  the  simulated  time 
trace  received  by  the  300th  receiver  matches  accurately 
with  the  experimentally  measured  one,  the  simulated  time 
trace  received  by  the  200th  receiver  is  substantially  differ¬ 
ent  compared  with  the  experimentally  measured  one.  To 
minimize  the  effects  of  model  mismatch,  we  replaced  these 
unreliable  measurements  with  computer-simulated  water 
bath  data,  as  described  in  Section  IV-C.  We  designated 
the  time  traces  received  by  the  512  receivers  located  on 
the  opposite  side  of  the  emitter  as  the  reliable  measure¬ 
ments  for  each  data  acquisition.  The  zeroth  data  acquisi¬ 
tion  of  the  combined  data  is  displayed  in  Fig.  18(b). 


(a)  (b) 

Fig.  19.  Images  reconstructed  from  the  experimentally  measured  phan¬ 
tom  data  by  use  of  (a)  the  bent-ray  model-based  sound  speed  recon¬ 
struction  method  and  (b)  the  WISE  method  with  a  TV  penalty  (/3TV 
=  1.0  x  102)  after  the  200th  iteration.  The  grayscale  window  is  [1.49, 
1.57]  mm/ps. 

Subsequently,  we  extracted  the  TOF  by  use  of  the  thresh¬ 
olding  method  with  a  thresholding  value  of  20%  of  the 
peak  value  of  each  time  trace.  The  bent-ray  reconstruction 
algorithm  was  applied  for  image  reconstruction  with  a 
measured  background  sound  speed  of  1.513  mm/ps.  The 
resulting  image  is  shown  in  Fig.  19(a)  and  has  a  pixel  size 
of  1  mm.  Finally,  the  image  was  smoothed  by  convolving 
it  with  a  2-D  Gaussian  kernel  with  a  standard  deviation 
of  2  mm. 


E.  Estimation  of  Initial  Guess 


F.  Image  Reconstruction 


The  initial  guess  for  the  WISE  method  was  obtained 
by  use  of  the  bent-ray  reconstruction  method  described  in 
Appendix  C.  We  first  filtered  each  time  trace  of  the  raw 
data  by  a  band-pass  Butterworth  filter  (0.5-2. 5  MHz). 


Fig.  18.  Zeroth  acquisition  of  (a)  the  experimentally  measured  raw  data 
and  (b)  the  combined  data,  and  time  traces  at  the  zeroth  acquisition 
received  by  (c)  the  300th  receiver,  and  (d)  the  200th  receiver.  The  gray¬ 
scale  window  for  (a)  and  (b)  is  [—45,  0]  dB. 


We  applied  the  WISE  method  with  a  TV  penalty  to  the 
combined  data  set.  The  second-order  wave  equation  solver 
was  employed  with  a  calculation  domain  of  dimensions 
512.0  x  512.0  mm2.  The  calculation  domain  was  sampled 
on  a  2560  x  2560  Cartesian  grid  with  a  grid  spacing  of 
0.2  mm.  On  a  platform  consisting  of  dual  quad-core  CPUs 
with  a  3.30  GHz  clock  speed,  64  GB  RAM,  and  a  single 
NVIDIA  Tesla  K20  GPU,  each  numerical  solver  run  took 
40  s  to  calculate  the  pressure  data  for  2112  time  samples. 
Knowing  the  size  of  the  phantom,  we  set  the  reconstruc¬ 
tion  region  to  be  within  a  circle  of  diameter  128  mm  (i.e., 
only  the  sound  speed  values  of  pixels  within  the  circle 
were  updated  during  the  iterative  image  reconstruction). 
We  swept  the  value  of  /2TV  over  a  wide  range  to  investigate 
its  impact  on  the  reconstructed  images. 

G.  Images  Reconstructed  From  Experimental  Data 

As  shown  in  Fig.  19,  the  spatial  resolution  of  the  image 
reconstructed  by  use  of  the  WISE  method  with  a  TV  pen¬ 
alty  is  significantly  higher  than  that  reconstructed  by  use 
of  the  bent-ray  model-based  method.  In  particular,  the 
structures  labeled  A  and  B  possess  clearly  defined  bound¬ 
aries.  This  observation  is  further  confirmed  by  the  profiles 
of  the  two  images  shown  in  Fig.  20.  In  addition,  the  struc¬ 
ture  labeled  C  in  Fig.  19(b)  is  almost  indistinguishable 
in  the  image  reconstructed  by  use  of  the  bent-ray  model- 
based  method  [Fig.  19(a)].  The  improved  spatial  resolu- 
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tion  is  expected  because  the  WISE  method  takes  into  ac¬ 
count  high-order  acoustic  diffraction,  which  is  ignored  by 
the  bent-ray  method  [23].  Though  not  shown  here,  for  the 
bent-ray  method,  we  investigated  multiple  TOF  pickers 
[25]  and  systematically  tuned  the  regularization  parame¬ 
ter.  As  such,  it  is  likely  that  Fig.  19(a)  represents  a  nearly 
optimal  bent-ray  image  in  terms  of  the  resolution.  This 
resolution  also  appears  to  be  similar  to  previous  experi¬ 
mental  results  reported  in  the  literature  [26]. 

The  convergence  properties  of  the  WISE  method  with 
a  TV  penalty  with  experimental  data  were  consistent  with 
those  observed  in  the  computer  simulation  studies.  Images 
reconstructed  by  use  of  10,  50,  and  300  algorithm  itera¬ 
tions  are  displayed  in  Fig.  21.  The  image  reconstructed  by 
use  of  10  iterations  contains  radial  streak  artifacts  that 
are  similar  in  nature  to  those  observed  in  the  computer 
simulation  studies.  These  artifacts  were  mitigated  after 
more  iterations.  The  image  reconstructed  after  300  itera¬ 
tions  [Fig.  21(d)]  appears  to  be  similar  to  that  after  200 
iterations  [Fig.  19(b)],  suggesting  that  the  WISE  method 
with  a  TV  penalty  is  close  to  convergence  after  about  200 
iterations.  The  time  required  to  complete  200  iterations 
was  approximately  14  hours.  The  estimated  time  it  would 
take  for  the  sequential  waveform  inversion  method  to  pro¬ 
duce  a  comparable  image  is  approximately  1  month,  as¬ 
suming  the  same  number  of  iterations  is  required  as  in  the 
computer  simulation  studies  (i.e.,  40). 

Despite  the  nonlinearity  of  the  WISE  method,  the  im¬ 
pact  of  the  TV  penalty  appears  to  be  similar  to  that  ob¬ 
served  in  other  imaging  applications  [54],  [58]  (Fig.  22). 
Though  not  shown  here,  the  impact  of  the  quadratic 
penalty  is  also  similar.  As  expected,  a  larger  value  of  (3 
reduced  the  noise  level  at  the  expense  of  spatial  image 
resolution.  These  results  suggest  a  predictable  impact  of 
the  penalties  on  the  images  reconstructed  by  use  of  the 
WISE  method. 


VII.  Summary 

It  is  known  that  waveform  inversion-based  reconstruc¬ 
tion  methods  can  produce  sound  speed  images  that  pos¬ 
sess  improved  spatial  resolution  properties  over  those 
produced  by  ray-based  methods.  However,  waveform  in¬ 
version  methods  are  computationally  demanding  and  have 
not  been  applied  widely  in  USCT  breast  imaging.  In  this 
work,  based  on  the  time-domain  wave  equation  and  mo¬ 
tivated  by  recent  mathematical  results  in  the  geophys¬ 
ics  literature,  the  WISE  method  was  developed  that  cir¬ 
cumvents  the  large  computational  burden  of  conventional 
waveform  inversion  methods.  This  method  encodes  the 
measurement  data  using  a  random  encoding  vector  and 
determines  an  estimate  of  the  sound  speed  distribution 
by  solving  a  stochastic  optimization  problem  by  use  of 
a  stochastic  gradient  descent  algorithm.  With  our  cur¬ 
rent  GPU-based  implementation,  the  computation  time 
was  reduced  from  weeks  to  hours.  The  WISE  method  was 
systematically  investigated  in  computer  simulation  and 


y  [mm] 

(a) 


Fig.  20.  Profiles  at  (a)  x  —  —24.0  mm  and  (b)  x  —  10.0  mm  of  the  recon¬ 
structed  images  by  use  of  the  bent-ray  model-based  sound  speed  recon¬ 
struction  method  (light  solid)  and  the  WISE  method  with  a  TV  penalty 
with  /3TV  =  1.0  x  102  (dark  dashed)  from  experimentally  measured  data. 

experimental  studies  involving  a  breast  phantom.  The 
results  suggest  that  the  method  holds  value  for  USCT 
breast  imaging  applications  in  a  practical  setting. 

Many  opportunities  remain  to  further  improve  the  per¬ 
formance  of  the  WISE  method.  As  shown  in  Fig.  19,  im¬ 
ages  reconstructed  by  use  of  the  WISE  method  can  con¬ 
tain  certain  artifacts  that  are  not  present  in  the  image 
reconstructed  by  use  of  the  bent-ray  method.  An  example 
of  such  an  artifact  is  the  dark  horizontal  streak  below 
the  structure  C.  Because  of  the  nonlinearity  of  the  im¬ 
age  reconstruction  problem,  it  is  challenging  to  determine 
whether  these  artifacts  are  caused  by  imaging  model  er¬ 
rors  or  by  the  optimization  algorithm,  which  might  have 
arrived  at  a  local  minimum  of  the  cost  function.  A  more 
accurate  imaging  model  can  be  developed  to  account  for 
out-of-plane  scattering,  the  transducer  finite  aperture  size 
effect,  acoustic  absorption,  as  well  as  other  physical  fac¬ 
tors.  Also,  the  stochastic  gradient  descent  algorithm  is 
one  of  the  most  basic  stochastic  optimization  algorithms. 
Numerous  emerging  optimization  algorithms  can  be  em¬ 
ployed  [33],  [42]  to  improve  the  convergence  rate.  In  ad- 
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(C)  (d) 

Fig.  21.  (a)  The  initial  guess  of  the  sound  speed  map  and  the  images 
reconstructed  by  use  of  the  WISE  method  with  a  TV  penalty  with  (/dTV 
=  1.0  x  102)  after  (b)  the  10th,  (b)  the  50th,  and  (d)  the  300th  iteration, 
from  the  experimentally  measured  phantom  data.  The  grayscale  window 
is  [1.49,  1.57]  mm/fis. 


dition,  there  remains  a  great  need  to  compare  the  WISE 
method  with  other  existing  sound  speed  reconstruction 
algorithms  [19],  [40]. 

There  remains  a  need  to  conduct  additional  investiga¬ 
tions  of  the  numerical  properties  of  the  WISE  method. 
Currently,  a  systematic  comparison  of  the  statistical  prop¬ 
erties  of  the  WISE  and  the  sequential  waveform  inversion 
method  is  prohibited  by  the  excessively  long  computa¬ 
tion  times  required  by  the  latter  method.  This  comparison 


will  be  interesting  when  a  more  efficient  wave  equation 
solver  is  available.  Given  the  fact  that  waveform  inver¬ 
sion  is  nonlinear  and  sensitive  to  the  initial  guess,  it  be¬ 
comes  important  to  investigate  how  to  obtain  an  accurate 
initial  guess.  We  also  observed  that  the  performance  of 
the  WISE  method  is  sensitive  to  how  strong  the  medium 
heterogeneities  are  and  the  profile  of  the  excitation  pulse. 
An  investigation  of  the  impact  of  the  excitation  pulse  the 
numerical  properties  of  the  image  reconstruction  may  help 
optimize  hardware  design.  In  addition,  quantifying  the 
statistics  of  the  reconstructed  images  will  allow  applica¬ 
tion  of  task-based  measures  of  image  quality  to  be  applied 
to  guide  system  optimization  studies. 


Appendix  A 

Continuous-to-Discrete  USCT  Imaging  Model 


In  practice,  each  data  function  gm(r,t)  is  spatially  and 

temporally  sampled  to  form  a  data  vector  gm  £  RNre°L, 
where  A1-60  and  L  denote  the  number  of  receivers  and  the 
number  of  time  samples,  respectively.  We  will  assume  that 
N*ec  and  L  do  not  vary  with  excitation  pulse.  Let  [gm]nvecL+l 
denotes  the  ( nrecL  +  l)th  element  of  gm.  When  the  receiv¬ 
ers  are  point-like,  gm  is  defined  as 

[gj„n+!  =  9m(r(m,niec)’1  A‘),  (30) 


where  the  indices  nrec  and  l  specify  the  receiver  location 
and  temporal  sample,  respectively,  and  A1  is  the  temporal 
sampling  interval.  The  vector  r(m,  nrec)  £  Q,m  denotes  the 
location  of  the  nrecth  receiver  at  the  rath  data  acquisition. 

A  C-D  imaging  model  for  USCT  describes  the  mapping 
of  c(r)  to  the  data  vector  gm  and  can  be  expressed  as 


fern] 


recL+l 


=  MmHcsm(r,  t) 


r=r(m,nTec),t=lAt 


for 


nrec  =  0,l,---,Arec  -1 

l  =  0, !,-••,£  —  1. 


(31) 


(a)  (b) 

Fig.  22.  Images  reconstructed  by  use  of  the  WISE  method  with  a  TV 
penalty  with  (a)  /7rv  m  5.0  x  101,  and  (b)  /3TV  =  5.0  x  102,  from  the 
experimentally  measured  phantom  data.  The  grayscale  window  is  [1.49, 
1.57]  mm/|is. 


Note  that  the  acousto-electrical  impulse  response  [59]  of 
the  receivers  can  be  incorporated  into  the  C-D  imaging 
model  by  temporally  convolving  sm(r,t)  in  (1)  with  the 
receivers’  acousto-electrical  impulse  response  if  we  assume 
all  receiving  transducers  share  an  identical  acousto-electri¬ 
cal  impulse  response. 

Appendix  B 

Frechet  Derivative  of  Data  Fidelity  Term 

Consider  the  integrated  squared-error  data  misfit  func¬ 
tion  [22],  [23], 
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where  gm(T,t)  and  gm(r,t)  denote  the  measured  data  func¬ 
tion  and  the  predicted  data  function  computed  by  use  of 
(3)  with  the  current  estimate  of  c(r). 

Both  the  sequential  and  WISE  reconstruction  method 
described  in  Section  III  require  knowledge  of  the  Frechet 
derivatives  of  JFcc(c)  and  7Zcc(c)  with  respect  to  c,  de¬ 
noted  by  V  cfFcc  and  Vc7^cc,  respectively.  The  calcula¬ 
tion  of  Vc7^cc  can  be  readily  accomplished  for  quadratic 
smoothness  penalties  [54] ,  [60] .  For  the  integrated  squared 
error  data  misfit  function  given  in  (32),  VCTCC  can  be 
computed  via  an  adjoint  state  method  as  [28],  [61],  [62] 

Vc^CC  =  FX  dtq™(r’T  ~  *)  (33) 

where  #m(r,t)  G  L2(M3  x  [0,oo))  is  the  solution  to  the  ad¬ 
joint  wave  equation.  The  adjoint  wave  equation  is  defined 
as 

1  d 2 

V  2qm( r,  t)  -  —  -2~  qm{ r,  t)  =  -rm(r,  t) ,  (34) 

c  (r)  a  t 

where  rjr,t)  =  gm(r,T  -  t)  -  gJr,T  -  t).  The  adjoint 
wave  equation  is  nearly  identical  in  form  to  the  wave  equa¬ 
tion  in  (1)  except  for  the  different  source  term  on  the 
right-hand  side,  suggesting  the  same  numerical  approach 
can  be  employed  to  solve  both  equations.  Because  one 
needs  to  solve  (1)  and  (34)  M  times  to  calculate  Vc^cc,  it 
is  generally  true  that  the  sequential  waveform  inversion  is 
computationally  demanding  even  for  a  2-D  geometry  [63] . 


Appendix  C 

Bent-Ray  Model-Based  Sound  Speed 
Reconstruction 

We  developed  an  iterative  image  reconstruction  algo¬ 
rithm  based  on  a  bent-ray  imaging  model.  The  bent-ray 
imaging  model  assumes  that  an  acoustic  pulse  travels 
along  a  ray  path  that  connects  the  emitter  and  the  re¬ 
ceiver  and  accounts  for  the  refraction  of  rays,  also  known 
as  ray-bending,  through  an  acoustically  inhomogeneous 
medium.  For  each  pair  of  receiver  and  emitter,  the  travel 
time,  as  well  as  the  ray  path,  is  determined  by  the  me¬ 
dium’s  sound  speed  distribution.  Given  the  travel  times 
for  a  collection  of  emitter  and  receiver  pairs  distributed 
around  the  object,  the  medium  sound  speed  distribution 
can  be  iteratively  reconstructed.  This  bent-ray  model- 
based  sound  speed  reconstruction  (BRSR)  method  has 
been  employed  in  the  USCT  literature  [26],  [64],  [65]. 

To  perform  the  BRSR,  we  extracted  a  TOF  data  vector 
from  the  measured  pressure  data.  Denoting  the  TOF  data 

vector  by  T  G  MMNre°,  each  element  of  T  represented  the 
TOF  from  each  emitter  and  receiver  pair.  The  extraction 
of  the  TOF  was  conducted  in  two  steps.  First,  we  esti¬ 
mated  the  difference  between  the  TOF  when  the  object 
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was  present  and  the  TOF  when  the  object  was  absent  by 
use  of  a  thresholding  method  [25],  [66].  In  particular,  20% 
of  the  peak  value  of  each  time  trace  was  employed  as  the 
thresholding  value.  Second,  a  TOF  offset  was  added  to  the 
estimated  difference  TOF  for  each  emitter  and  receiver 
pair  to  obtain  the  absolute  TOF,  where  the  TOF  offset 
was  calculated  according  to  the  scanning  geometry  and 
the  known  background  SOS. 

Having  the  TOF  vector  T,  we  reconstructed  the  sound 
speed  by  solving  the  following  optimization  problem: 

s  =  argmin||T  —  Kss||2  +  /37Z(s) ,  (35) 

s 

where  s  denotes  the  slowness  (the  reciprocal  of  the  SOS) 
vector,  and  Ks  denotes  the  system  matrix  that  maps  the 
slowness  distribution  to  the  TOF  data.  The  superscript  s 
indicates  the  dependence  of  Ks  on  the  slowness  map.  At 
each  iteration,  using  the  current  estimate  of  the  SOS,  a 
ray-tracing  method  [50]  was  employed  to  construct  the 
system  matrix  Ks.  Explicitly  storing  the  system  matrix  in 
a  sparse  representation,  we  utilized  the  limited  BFGS 
method  [51]  to  solve  the  optimization  problem  given  in 
(35).  The  estimated  slowness  was  then  converted  to  the 
sound  speed  by  taking  the  reciprocal  of  the  s  element.  We 
refer  the  readers  to  [26],  [64]— [67]  for  more  details  about 
the  BRSR  method. 
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ABSTRACT 

In  this  work,  we  investigate  a  novel  reconstruction  method  for  laser-induced  ultrasound  computed  tomography 
(USCT)  breast  imaging  that  circumvents  limitations  of  existing  methods  that  rely  on  ray-tracing.  There  is 
currently  great  interest  in  developing  hybrid  imaging  systems  that  combine  optoacoustic  tomography  (OAT)  and 
USCT.  There  are  two  primary  motivations  for  this:  (1)  the  speed-of-sound  (SOS)  distribution  reconstructed  by 
USCT  can  provide  complementary  diagnostic  information;  and  (2)  the  reconstructed  SOS  distribution  can  be 
incorporated  in  the  OAT  reconstruction  algorithm  to  improve  OAT  image  quality.  However,  image  reconstruction 
in  USCT  remains  challenging.  The  majority  of  existing  approaches  for  USCT  breast  imaging  involve  ray¬ 
tracing  to  establish  the  imaging  operator.  This  process  is  cumbersome  and  can  lead  to  inaccuracies  in  the 
reconstructed  SOS  images  in  the  presence  of  multiple  ray-paths  and/or  shadow  zones.  To  circumvent  these 
problems,  we  implemented  a  partial  differential  equation-based  Eulerian  approach  to  USCT  that  was  proposed 
in  the  mathematics  literature  but  never  investigated  for  medical  imaging  applications.  This  method  operates 
by  directly  inverting  the  Eikonal  equation  without  ray-tracing.  A  numerical  implementation  of  this  method  was 
developed  and  compared  to  existing  reconstruction  methods  for  USCT  breast  imaging.  We  demonstrated  the 
ability  of  the  new  method  to  reconstruct  SOS  maps  from  TOF  data  obtained  by  a  hybrid  OAT/USCT  imager 
built  by  our  team. 

Keywords:  ultrasound  tomography,  optoacoustic  tomography,  photoacoustic  tomography,  breast  cancer  imag¬ 
ing 


1.  INTRODUCTION 

Transmission  ultrasound  computed  tomography  (USCT)  is  an  emerging  imaging  modality  with  many  biomedical 
applications.  USCT  can  be  employed  to  retrieve  anatomical  information  of  tissues  e.  g.  speed  of  sound,  acoustical 
impedance  and  reflectivity.  The  effectiveness  of  USCT  in  tumor  detection  has  been  discussed  in  recent  studies.1-3 
It  is  known  that  cancerous  tissues  have  higher  SOS  values  compared  to  the  benign  fatty  masses  and  healthy 
breast  tissues.  A  clinical  ultrasound  ring  array  scanner  for  breast  cancer  diagnosis  (Computed  Ultrasound  Risk 
Evaluation  (CURE))  has  been  proposed.1,2  This  system  consists  of  256  transducer  elements  distributed  on 
a  ring  with  a  20  cm  diameter.  Another  prototype  has  also  been  developed  by  Soft  Vue  and  consists  of  2048 
transducers.4  The  system  is  capable  of  reconstructing  a  series  of  2D  slices  of  the  SOS,  acoustic  attenuation,  and 
reflectivity  distributions.  Techniscan  (Salt  Lake  City  UT)  introduced  a  commercial  USCT  system  that  employs 
three  transducer  probes  placed  around  the  breast.  The  transducer  system  is  mechanically  rotated  to  reconstruct 
2D  slices  and  subsequently  vertically  scanned  to  capture  multiple  slices  to  obtain  3D  images.5  Finally,  an 
ultrasound  imaging  module  capable  of  generating  three-dimensional  SOS  distributions  has  been  investigated  by 
researchers  at  the  Karlsruhe  Institute  of  Technology  (KIT),  Germany.6 

Biomedical  applications  of  USCT  commonly  employ  geometrical  acoustics  models  and  require  time-of-flight 
(TOF)  measurements.  TOF  data  that  have  been  recorded  for  many  source- receiver  pairs  can  be  employed  for 
reconstruction  of  the  SOS  distribution.  The  reconstruction  of  the  speed  of  sound  distribution  is  conventionally 
performed  by  using  ray-tracing  (RT)  methods.1,2, 7,8  To  account  for  the  curvature  of  the  ray  paths,  the  rays  are 
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traced  along  the  negative  gradient  of  the  TOF  distribution.8  Ray-tracing  can  become  cumbersome,  especially 
for  three-dimensional  USCT.  Moreover,  unlike  X-ray  computed  tomography,  the  heterogeneous  SOS  distribution 
results  in  uneven  ray-distributions,  which  makes  the  inverse  problem  ill-conditioned.  In  this  work,  we  will 
investigate  a  different  approach  for  SOS  image  reconstruction.  This  method,  the  adjoint  state  (AS)  method,  has 
previously  been  employed  for  seismic  tomography.9, 10  We  are  investigating  the  AS  method  for  USCT  for  the 
first  time  for  biomedical  applications. 

2.  DESCRIPTION  OF  NUMERICAL  STUDIES 

To  perform  USCT  using  RT,  we  have  employed  the  the  geometrical  ray  theory  for  sound  waves.  This  approxi¬ 
mation  results  in  a  non-linear  model,  the  eikonal  equation,  to  relate  the  measured  TOF  values  to  the  speed  of 
sound  distribution  as 


|VT(r)| 


1 

c(r)' 


(1) 


In  Eq.  (1),  VT  is  the  gradient  of  the  TOF,  T,  and  c  is  the  SOS  distribution,  both  of  which  are  a  function  of 
position  as  denoted  by  the  position  vector  r  £  M2.  Currently,  the  bent-ray  reconstruction  is  a  widely  employed 
reconstruction  technique  for  USCT  because  it  incorporates  refraction  during  sound  wave  propagation.  The 
eikonal  equation  is  solved  numerically  by  finite  difference  methods11  to  obtain  a  TOF  map  for  a  given  source 
corresponding  to  a  certain  speed  of  sound  map  c(r). 


2.1  Ray-tracing  reconstruction  method 

In  RT  methods,  the  TOF  is  calculated  as  the  line-integral  over  the  slowness  distribution  over  the  ray-path 
connecting  the  source  and  the  receiver  location: 

T{  r)  =  /  -b  (2) 

Jr(c)  c(r) 

The  dependence  of  the  ray-path,  T(c),  on  the  SOS  distribution  makes  it  a  non-linear  problem.  The  discretized 
imaging  model  is  given  by 


T  =  H(c)I  (3) 

where  T  is  a  vector  of  TOF  measurments,  c  is  a  finite-dimensional  representation  of  the  SOS,  and  H(c)  is  the 
system  matrix.  To  formulate  H(c),  we  implemented  a  RT  method.  Weights  were  assigned  to  pixels  in  the  discrete 
SOS  map  based  on  the  number  of  times  each  pixel  was  intersected  by  the  rays.  This  weight  matrix  constitutes 
H(c). 

To  estimate  the  SOS  distribution  from  the  measured  TOF  data,  we  solved  the  following  optimization  problem: 

c  =  argmin  ||  T  —  T*  ||2  +vg( c),  (4) 

c 

where  c  denotes  the  sought-after  estimate  of  the  the  SOS  distribution,  T*  is  the  measured  TOF  data  from 
all  source-transducer  pairs,  g( c)  is  a  penalty  function,  v  is  a  regularization  parameter,  and  T  is  the  computed 
TOF  found  by  solving  (1).  To  minimize  Eq.  (4)  we  used  the  Limited  BFGS  method.12  In  solving  the  nonlinear 
optimization  problem,  we  evaluated  the  gradient  of  Eq.  (4)  as 

Vc  =  2H(c)t  [H(c)c  -  T*]  +  vVg( c).  (5) 

It  should  be  noted  from  the  above  equation  that  the  first  term  is  a  linear  approximation  of  the  true  non-linear 
gradient  of  the  objective  function.  The  above  linearized  gradient  is  primarily  used  in  USCT  for  biomedical 
applications. 
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2.2  Adjoint-state-based  reconstruction  method 

We  implemented  a  previously  proposed  algorithm  for  USCT  reconstructed  based  on  the  adjoint-state  method.9 
The  mismatch  energy  functional  between  the  measured  and  simulated  data  is  defined  as9 

E[c(r)}  =  ±JjT(r)-T*(v)\2<m,  (6) 

where  T*|s  is  the  measured  TOF  and  T\s  is  computed  by  solving  Eq.  (1).  The  quantity  in  Eq.  (6)  —  the 
energy  functional  —  measures  the  L2-difference  between  the  the  solution  of  the  eikonal  equation,  T,  and  the 
experimental  measurement,  T*,  on  the  measurement  surface  S.  Using  the  adjoint-state  method  for  a  small 
perturbation  ec  to  c,  the  gradient  of  the  energy  is  defined: 

(7) 

Here,  V  is  volume  enclosed  the  measurement  surface  S  and  A(r)  is  the  adjoint  function  to  T(r)  that  satisfies  the 
following  adjoint  equation: 


V  •  [A(r)VT(r)]  =  0 


with  the  boundary  condition, 


[n-VT(r)]A(r)|s  =  [T*(r)-T(r)]s. 


(8) 

(9) 


Here  n  is  the  unit  outward  normal  of  the  surface  S.  To  minimize  the  energy  using  the  method  of  gradient 
descent,  a  perturbation  c(r)  =  —  A(r)/c3(r)  is  defined.  This  leads  to 

SE  =  —e  [  c2(r)dO  <  0,  (10) 

Jv 

where  V  denotes  the  region  interior  to  S.  By  solving  Eqs.  (8)  and  (9),  the  update,  c(r),  to  the  SOS  distribution 
can  be  computed.  Specifically,  the  SOS  distribution  is  updated  at  each  step  as: 

Ck+ 1  =Ck+6kCk  (11) 

until  a  convergence  criterion  is  reached.  The  following  two  conditions  are  required  of  the  SOS  distribution: 
(i)  ck\s  =  0  and  (ii)  ck+1  is  smooth.  To  fulfill  (ii),  a  regularization  term  similar  to  the  one  used  in  Eq.  (4) 
was  included.  The  filtering  scheme  defined  in  Lueng’s  work9  was  also  implemeted.  The  step  size,  ek ,  can  be 
determined  by  using  the  Armijo- Golstein  rule  or  by  simply  setting  ek  =  e.  The  update  scheme  described  in 
Eq.  (11)  takes  a  large  number  of  iterations  to  converge.  Therefore,  we  used  the  limited-memory  Broydon- 
Fletcher  Goldfarb-Shanno  (L-BFGS)  method  to  solve  for  this  nonlinear  optimization  problem.  We  solved  the 
adjoint-state  equation  (8)  using  the  fast-sweeping  method9  with  the  boundary  conditions  defined  in  Eq.  (9). 


2.3  Experimental  Setup 

The  experimental  setup  consists  of  a  single  laser  ultrasound  (LU)  source  and  a  64  transducer  elements  arranged 
in  an  arc.  The  array  aperture  spans  a  152  degree  arc  with  a  radius  of  65  mm.  The  imaging  module  is  mounted 
and  centered  on  a  rotational  stage  operated  by  a  stepper  motor,  which  is  used  to  obtain  TOF  measurement  for 
150  views.  The  distance  between  the  central  element  of  the  arc  array  and  the  LU  source  is  130  mm.  Optoacoustic 
imaging  was  concurrently  performed.  More  detail  about  the  LU  sources  and  the  transducer  array  can  be  found 
elsewhere.13, 14 
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Figure  1.  (a)  Speed  of  sound  distribution  for  the  phantom;  reconstructed  SOS  distribution  for  the  RT  method  (b)  and  for 
the  AS  method  (c). 


3.  RESULTS  AND  DISCUSSION 

Computer-simulation  studies  were  conducted  to  compare  the  RT  method  with  the  AS  method.  In  this  study,  the 
objective  function  did  not  include  a  penalty  and  least  squares  estimates  of  the  SOS  distributions  were  computed. 
The  two-dimensional  phantom  depicting  the  SOS  distribution  is  shown  in  Figure  1(a).  The  phantom  consists 
of  three  discs  of  4.72  mm  diameter  with  constant  SOS  of  1.48,  1.6  and  1.7  mm/gs,  respectively.  The  pressure 
data  were  generated  using  the  k-Wave  software  package15  with  a  geometry  and  acoustic  properties  consistent 
with  our  experimental  system  design.  Gaussian  white  noise  was  added  to  the  calculated  pressure  signal  to  obtain 
experimentally-relevant  SNRs.  Figures  1(b)  and  1(c)  show  the  reconstructed  SOS  distribution  for  the  RT  method 
and  the  AS  method,  respectively.  In  the  case  of  the  RT  method,  streak  artifacts  are  very  visible  as  compared  to 
the  AS  method.  All  three  structures  can  be  seen  in  the  AS  reconstruction  of  the  SOS  distribution.  The  results 
show  that  the  AS  method  can  be  successfully  used  to  perform  USCT  for  biomedical  applications. 

To  check  the  accuracy  of  the  reconstructed  SOS  distribution,  we  selected  2  mm  x  2mm  regions  at  different 
locations  and  calculated  the  avaeraged  SOS  and  standard  deviation  in  those  region.  The  location  of  the  five 
regions  is  shown  in  Fig. 2 (a).  The  bar  plot  in  Fig. 2(b)  shows  the  averaged  SOS  values  for  both  the  RT  and  AS 
method.  It  can  be  seen  that  AS  gives  accurate  SOS  values  for  the  selected  regions.  The  maximum  standard 
deviation  was  0.0226  mm/gs  for  region  ”A”. 

Finally,  we  studied  the  use  of  the  AS  method  to  perform  SOS  image  reconstruction  for  the  experimentally 
measured  TOF  from  our  LU  system.  The  experimental  phantom  consists  of  three  tubes  each  with  a  4.72  mm 
internal  diameter.  Tubes  were  filled  with  water  at  different  salt  concentrations  to  induce  different  SOS  values 
and  CUSO4  was  added  to  provide  optical  contrast.  For  this  case,  concurrent  optoacoustic  (OA)  and  ultrasonic 
data  acquisition  was  performed.  The  OA  reconstruction  of  the  phantom  is  shown  in  Fig.  3(a)  and  the  SOS 
reconstruction  is  shown  in  Fig.  3(b).  The  SOS  image  was  found  via  the  AS  reconstruction  method.  Once  again, 
all  three  discs  are  visible  in  the  reconstructed  SOS  distribution. 


I  Phantom 
]AS  method 
IRT  method 


Figure  2.  (a)  Speed  of  sound  distribution  for  the  phantom  with  the  marked  location  of  five  regions;  (b)  bar  plot  for  the 
averaged  SOS  values  in  five  regions  for  phantom,  AS  method,  and  RT  method. 
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1.62 


Figure  3.  (a)  Optoacoustic  image  of  three  tubes;  (b)  reconstructed  SOS  distribution  using  the  AS  method. 

4.  SUMMARY 

The  adjoint  state  method  has  been  implemented  for  biomedical  applications  of  USCT.  Images  reconstructed  from 
both  simulation  studies  and  measured  TOF  data  were  presented.  Ray  tracing  becomes  much  more  cumbersome 
for  three-dimensional  USCT.  Consequently,  the  adjoint  state  method  holds  great  promise  for  that  application. 
Further  numerical  studies  will  also  be  performed  to  quantify  resolution  and  noise  propagation  in  the  adjoint  state 
method. 
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ABSTRACT 


In  this  work,  we  investigate  a  novel  reconstruction  method  for  laser-induced 
ultrasound  tomography  (UST)  breast  imaging  that  circumvents  limitations  of 
existing  methods  that  rely  on  ray-tracing.  There  is  currently  great  interest  in 
developing  hybrid  imaging  systems  that  combine  optoacoustic  tomography 
(OAT)  and  UST.  There  are  two  primary  motivations  for  this:  (1)  the  speed- 
of-sound  (SOS)  distribution  reconstructed  by  UST  can  provide 
complementary  diagnostic  information;  and  (2)  the  reconstructed  SOS 
distribution  can  be  incorporated  in  the  OAT  reconstruction  algorithm  to 
improve  OAT  image  quality.  However,  image  reconstruction  in  UST 
remains  challenging.  The  majority  of  existing  approaches  for  UST  breast 
imaging  involve  ray-tracing  to  establish  the  imaging  operator.  This  process 
is  cumbersome  and  can  lead  to  severe  inaccuracies  in  the  reconstructed  SOS 
images  in  the  presence  of  multiple  ray-paths  and/or  shadow  zones. 

To  circumvent  these  problems,  we  implemented  a  partial  differential 
equation-based  Eulerian  approach  to  UST  that  was  proposed  in  the 
mathematics  literature  but  never  investigated  for  medical  imaging 
applications.  This  method  operates  by  directly  inverting  the  Eikonal 
equation  without  ray-tracing.  A  numerical  implementation  of  this  method 
was  developed  and  systematically  compared  to  existing  reconstruction 
methods  for  UST  breast  imaging.  We  demonstrated  the  ability  of  the  new 
method  to  reconstruct  accurate  SOS  maps  from  TOF  data  obtained  by  a  3D 
hybrid  OAT/UST  imager  built  by  our  team. 
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ABSTRACT 

Ultrasound  computed  tomography  (USCT)  holds  great  promise  for  improving  the  detection  and  management  of 
breast  cancer.  Because  they  are  based  on  the  acoustic  wave  equation,  waveform  inversion-based  reconstruction 
methods  can  produce  images  that  possess  improved  spatial  resolution  properties  over  those  produced  by  ray- 
based  methods.  However,  waveform  inversion  methods  are  computationally  demanding  and  have  not  been  applied 
widely  in  USCT  breast  imaging.  A  computationally  efficient  numerical  wave  equation  solver  has  been  reported 
based  on  a  modified  Fresnel  propagation,  which  only  applies  to  USCT  systems  with  a  planar  incident  wave.  For 
breast  imaging  systems  with  a  spherical  incident  wave,  waveform  inversion-based  reconstruction  methods  remain 
computationally  challenging. 

In  this  work,  source  encoding  concepts  are  employed  to  develop  an  accelerated  USCT  reconstruction  method 
that  circumvents  the  large  computational  burden  of  conventional  waveform  inversion  methods.  This  method, 
referred  to  as  the  waveform  inversion  with  source  encoding  (WISE)  method,  encodes  the  measurement  data 
using  a  random  encoding  vector  and  determines  an  estimate  of  the  speed-of-sound  distribution  by  solving  a 
stochastic  optimization  problem  by  use  of  a  stochastic  gradient  descent  algorithm.  For  practical  applications,  a 
data- filling  strategy  is  proposed  to  mitigate  source  inferences  to  its  neighbor  receivers.  Computer-simulation  and 
experimental  phantom  studies  are  conducted  to  demonstrate  the  use  of  the  WISE  method.  Using  a  single  graphics 
processing  unit  card,  each  iteration  can  be  completed  within  25  seconds  for  a  128  x  128  mm2  reconstruction  region. 
The  results  suggest  that  the  WISE  method  maintains  the  high  spatial  resolution  of  waveform  inversion  methods 
while  significantly  reducing  the  computational  burden. 

1.  PURPOSE 

This  study  is  focused  on  the  image  reconstruction  of  breast  speed-of-sound  (SOS)  distribution  in  USCT.  The 
majority  of  USCT  image  reconstruction  methods  for  breast  imaging  investigated  to  date  have  been  based  on 
approximations  to  the  acoustic  wave  equation.1,2  A  relatively  popular  class  of  methods  is  based  on  geometrical 
acoustics.  They  are  commonly  referred  to  as  ‘ray-based’  methods.  Although  ray-based  methods  can  be  compu¬ 
tationally  efficient,  the  spatial  resoultion  of  the  images  they  produce  is  limited  due  to  the  fact  that  diffraction 
effects  are  not  modelled.3,4  This  is  undesirable  for  breast  imaging  applications,  in  which  the  ability  to  resolve 
fine  features,  e.g.,  tumor  spiculations,  is  important  for  distinguishing  healthy  from  diseased  tissues. 

USCT  reconstruction  methods  based  on  the  acoustic  wave  equation,  also  known  as  full-wave  inverse  scattering 
or  waveform  inversion  methods,  have  also  been  explored  for  a  variety  of  applications  including  medical  imaging.4-7 
Because  they  account  for  higher-order  diffraction  effects,  waveform  inversion  methods  can  produce  images  that 
possess  higher  spatial  resolution  properties  than  those  produced  by  ray-based  methods.4,5  However,  conventional 
waveform  inversion  methods  are  iterative  in  nature  and  require  the  wave  equation  to  be  solved  numerically  a  large 
number  of  times  at  each  iteration.  Consequently,  such  methods  can  be  extremely  computationally  burdensome. 
For  special  geometries,7  efficient  numerical  wave  equation  solvers  have  been  reported.  However,  apart  from 
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special  cases,  the  large  computational  burden  of  waveform  inversion  methods  has  hindered  their  widespread 
application. 

The  purpose  of  this  study  is  to  develop  an  algorithmically  accelerated  waveform  inversion  method  for  breast 
SOS  reconstruction.  Aided  by  a  graphics  processing  unit  (GPU)-accelerated  implementation,  the  developed 
method  will  maintain  the  high  spatial  resolution  of  standard  waveform  inversion  methods  with  a  signifcant 
reduction  in  computational  time. 


2.  METHODS 

A  conventional  waveform  inversion  method  seeks  the  solution  of 

1  M—l 

c  =  argmin  -  V  ||gm  -  Hcsm||2  +  f3H{ c),  (1) 

c  Z  z '  - 

m— 0 

where  c  is  the  sought-after  object  to  be  reconstructed,  i.e,  SOS  distribution,  gm  denotes  the  measured  data 
vector,  sm  denotes  the  (known)  source  vector,  Hc  denotes  a  numerical  wave  equation  solver  (NWES)  that 
maps  the  known  source  vector  to  the  measured  data  vector,  and  71(c)  and  (3  denote  the  penalty  term  and  the 
regularization  parameter  respectively.  The  superscript  in  Hc  indicates  the  dependence  of  Hc  on  c.  Note  that  one 
USCT  measurement  involves  firing  a  sequence  of  acoustic  pulses  in  turn  and  recording  the  data  corresponding  to 
every  pulse.  Each  pulse- firing  and  data  recording  process  will  be  indexed  by  m  for  m  =  0, 1,  •  •  •  ,  M  —  1.  Solving 
Eqn.  (1),  in  general,  requires  the  calculation  of  \  Vc||gm_ Hcsm||2,  where  Vc  denotes  the  gradient  operator 

with  respect  to  c.  The  gradient  in  each  summand  is  commonly  computed  by  use  of  an  adjoint  state  method,5 
which  requires  two  runs  of  the  NWES.  Repeating  the  gradient  calculation  for  all  sources  results  in  2 M  runs 
of  the  NWES  at  each  iteration.  This  computational  burden  largely  hinders  the  application  of  the  conventional 
waveform  inversion  methods  in  practice. 

In  this  study,  a  waveform  inversion  with  source  encoding  (WISE)  method  was  developed.  The  WISE  method 
employs  the  objective  function 


c  =  argminEw{I||gw  -  Hcsw||2}  +  (3K(c),  (2) 

where  Ew  denotes  the  expectation  operator  with  respect  to  the  random  source  encoding  vector  w  G  Mm,  and 
gw  and  sw  denote  the  w-encoded  data  and  source  vectors,  defined  as 


M-l  M  —  l 

gw  =  [w]mgm,  and  sw=  y][w]msm,  (3) 

m=0  m=0 

respectively.  Equation  (2)  was  solved  by  use  of  a  stochastic  gradient  descent  algorithm.8  Because  the  stochastic 
gradient  descent  algorithm  calculated  the  gradient  of  only  one  realization  of  the  random  variable  \  1 1  gw — MHcsw||2 
at  each  iteration,  the  required  number  of  NWES  runs  per  iteration  was  reduced  from  2 M  to  2.  Although  it,  in 
general,  requires  more  algorithm  iterations  to  average  out  the  randomness  in  the  realizations,  the  WISE  method, 
as  demonstrated  later,  can  greatly  reduce  the  overall  number  of  NWES  runs.  Both  computer-simulation  and 
experimental  phantom  studies  were  conducted  to  demonstrate  the  use  of  the  WISE  method  for  breast  SOS 
reconstruction. 


3.  RESULTS 

The  images  reconstructed  from  the  computer-simulated  noise- free  data  by  use  of  the  WISE  method  after  199 
iterations  and  sequential  waveform  inversion  method  after  43  iterations  are  shown  in  Fig.  l-(a)  and  (b).  As 
expected,4,9  both  images  are  more  accurate  and  possess  higher  spatial  resolution  than  the  one  reconstructed  by 
use  of  the  bent-ray  reconstruction  algorithm  displayed  in  Fig.  l-(c).  The  images  shown  in  Fig.  l-(a)  and  -(b) 
possess  similar  accuracies  as  measured  by  their  Euclidean  distances  from  the  SOS  phantom  vector  c,  namely 
0.07%  of  ||c||  for  the  former  and  0.08%  of  ||c||  for  the  latter.  However,  the  reconstruction  of  Fig.  l-(a)  required 


only  about  1.7%  of  the  computational  time  required  to  reconstruct  Fig.  1- (b) ,  namely,  1.4  hours  for  the  former 
and  81.4  hours  for  the  latter  respectively.  This  is  because  the  WISE  method  required  only  1018  NWES  runs, 
which  is  signficantly  less  than  the  58880  NWES  runs  required  by  the  sequential  waveform  inversion  method. 
With  a  similar  number  of  NWES  runs,  (e.g.,  1024),  one  can  only  complete  a  single  algorithm  iteration  by  use 
of  the  sequential  waveform  inversion  method.  The  corresponding  image,  shown  in  Fig.  1- (d) ,  lacks  quantitative 
accuracy  as  well  as  qualitative  value  for  identifying  features.  The  results  suggest  that  the  WISE  method  maintains 
the  advantages  of  the  sequential  waveform  inversion  method  while  significantly  reducing  the  computational  time. 


(a)  (b)  (c)  (d) 

Figure  1.  Images  reconstructed  by  use  of  (a)  the  WISE  method  after  the  199-th  iteration  (1,018  runs  of  NWES)  (b)  the 
sequential  waveform  inversion  algorithm  after  the  43-rd  iteration  (58,880  runs  of  the  NWES),  (c)  the  bent-ray  model- 
based  SOS  reconstruction  method,  and  (d)  the  sequential  waveform  inversion  algorithm  after  the  1-st  iteration  (1,024 
runs  of  the  NWES)  from  the  noise-free  non- attenuated  data.  The  grayscale  window  is  [1.46,  1.58]  mm//xs. 


The  images  reconstructed  from  the  experiment  ally- measured  data  are  shown  in  Fig.  2.  The  spatial  resolution 
of  the  image  reconstructed  by  use  of  the  WISE  method  is  significantly  higher  than  that  reconstructed  by  use 
of  the  bent-ray  model-based  method.  In  particular,  the  structures  labeled  ‘A’  and  CB’  possess  clearly- defined 
boundaries.  In  addition,  the  structure  labeled  ‘Cancer’  in  Fig.  2- (a)  is  almost  indistinguishable  in  the  image 
reconstructed  by  use  of  the  bent-ray  model-based  method  (see  Fig.  2-(b)).  The  improved  spatial  resolution  is 
expected  because  the  WISE  method  takes  into  account  the  high-order  diffraction  effects,  which  are  ignored  by 
the  bent-ray  method.4 


Fibroadenoma 


Figure  2.  (a)  Schematic  of  the  breast  phantom  employed  in  the  experimental  study.  Images  reconstructed  from  the 
experimentally  measured  phantom  data  by  use  of  (a)  the  bent-ray  model-based  SOS  reconstruction  method  and  (b)  the 
WISE  method  after  the  200-th  iteration.  The  grayscale  window  is  [1.49, 1.57]  mm//is. 


4.  NEW  OR  BREAKTHROUGH  WORK  TO  BE  PRESENTED 

Source  encoding  concepts  are  demonstrated  in  breast  USCT  experimental  studies  for  the  first  time.  Unlike 
previously  studied  waveform  inversion  methods  that  were  based  on  the  Helmholtz  equation,  the  WISE  method 


is  formulated  by  use  of  the  time-domain  acoustic  wave  equation.  A  GPU-accelerated  NWES  is  developed  that 
can  compute  1800  time  samples,  on  a  1024  x  1024  spatial  grid,  in  5  seconds.  In  addition,  a  data-filling  strategy 
is  proposed  to  mitigate  the  inference  of  the  source  with  its  neighboring  receivers  for  practical  applications. 

5.  CONCLUSION 

It  is  known  that  waveform  inversion-based  reconstruction  methods  can  produce  SOS  images  that  possess  im¬ 
proved  spatial  resolution  properties  over  those  produced  by  ray-based  methods.  However,  waveform  inversion 
methods  are  computationally  demanding  and  have  not  been  applied  widely  in  USCT  breast  imaging.  In  this 
work,  based  on  the  time-domain  wave  equation  and  motivated  by  recent  mathematical  results  in  the  geophysics 
literature,  the  WISE  method  was  developed  that  circumvents  the  large  computational  burden  of  conventional 
waveform  inversion  methods.  This  method  encodes  the  measurement  data  using  a  random  encoding  vector  and 
determines  an  estimate  of  the  speed-of-sound  distribution  by  solving  a  stochastic  optimization  problem  by  use 
of  a  stochastic  gradient  descent  algorithm.  With  our  current  GPU-based  implementation,  the  computation  time 
was  reduced  from  weeks  to  hours.  The  WISE  method  was  systematically  investigated  in  computer- simulation 
and  experimental  studies  involving  a  breast  phantom.  The  results  suggest  that  the  method  holds  value  for  USCT 
breast  imaging  applications  in  a  practical  setting. 

6.  DISCLOSURE 

This  work  is  original.  Parts  of  this  work  have  been  submitted  to  IEEE  Transactions  on  Ultrasounics,  Ferro- 
electrics  and  Frequency  Control  and  are  under  review. 
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ABSTRACT 

Iterative  image  reconstruction  algorithms  can  model  complicated  imaging  physics,  compensate  for  imperfect 
data  acquisition  systems,  and  exploit  prior  information  regarding  the  object.  Hence,  they  produce  higher  quality 
images  than  do  analytical  image  reconstruction  algorithms.  However,  three-dimensional  (3D)  iterative  image 
reconstruction  is  computationally  burdensome,  which  greatly  hinders  its  use  with  applications  requiring  a  large 
field-of-view  (FOV),  such  as  breast  imaging.  In  this  study,  an  improved  GPU-based  implementation  of  a  nu¬ 
merical  imaging  model  and  its  adjoint  have  been  developed  for  use  with  general  gradient-based  iterative  image 
reconstruction  algorithms.  Both  computer  simulations  and  experimental  studies  are  conducted  to  investigate  the 
efficiency  and  accuracy  of  the  proposed  implementation  for  optoacoustic  tomography  (OAT).  The  results  suggest 
that  the  proposed  implementation  is  more  than  five  times  faster  than  the  previous  implementation. 

Keywords:  Optoacoustic  tomography,  iterative  image  reconstruction,  GPU  acceleration 

1.  INTRODUCTION 

Optoacoustic  tomography  (OAT)  is  an  emerging  imaging  modality  with  many  biomedical  applications.1-3  The 
clinical  applications  of  OAT,  e.g.,  breast  imaging,  have  been  pursued  by  several  research  groups.4-6  OAT  imaging 
is  performed  by  illuminating  tissues  with  a  laser  pulse.  The  thermoacoustic  effect  then  creates  an  acoustic  pressure 
rise  that  can  be  measured  by  transducer  elements  placed  around  the  object.1,7  From  these  measured  pressure 
data,  we  can  reconstruct  an  estimate  of  the  absorbed  optical  energy  density  of  the  object.  Many  analytical  and 
iterative  image  reconstruction  algorithms  have  been  developed  for  this  purpose.8  Analytical  methods,  such  as 
filtered  backprojection  (FBP),8-10  have  the  advantages  of  being  simplistic  and  computationally  efficient;  however, 
due  to  model  errors,  measurement  noise,  data  inconsistencies,  and  limited  data,11  these  methods  often  produce 
lower  quality  images  than  those  obtained  by  iterative  methods.  Iterative  algorithms  have  also  been  used  to 
incorporate  a  more  accurate  imaging  physics  as  well  as  the  spatial  and  acousto-electric  impulse  responses  of 
ultrasonic  transducers.11-13 

Despite  the  many  advantages  of  iterative  methods,  few  studies  have  been  reported  involving  3D  iterative 
image  reconstruction  algorithms,12,14,15  primarily  because  of  their  computational  complexity.  Progress  has  been 
made  in  this  regard  by  the  use  of  graphic  processing  unit  (GPU)  accelerated  algorithms;  however,  even  in  this 
case,  it  takes  tens  of  minutes  to  perform  a  single  iteration  and  several  hours  to  obtain  an  image,  making  iterative 
algorithms  less  attractive  for  practical  applications.15  In  this  study,  we  present  an  accelerated  interpolation-based 
imaging  algorithm  that  combines  four  strategies:  (i)  the  use  of  a  new  imaging  model  that  eliminates  one  operator, 
(ii)  the  use  of  a  voxel-driven  adjoint  imaging  operator,  (iii)  the  acceleration  of  convolution  calculations  by 
CUDA’s  parallelized  fast-Fourier  transform  library,  and  (iv)  a  multi-GPU  implementation.  Image  reconstruction 
was  performed  by  solving  a  panelized  least  squares  optimization  problem  using  a  linear  conjugate  gradient 
algorithm.12, 16 
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2.  DESCRIPTION  OF  THE  IMAGING  MODEL 


2.1  Continuous-to-Discrete  (C-D)  Imaging  Model 

When  sampling  effects  are  considered,  an  OAT  system  for  a  point-like  transducer  is  properly  described  by  a 
continuous-to-discrete  (C-D)  imaging  model: 


U]  qK+k 


de(t)*tg(r8,t) 


q= 0,1,---  ,Q- 1 
fe=0,l,---  ,K-V 


(1) 


where  K  represents  the  total  number  of  temporal  samples,  Q  is  the  number  of  transducer  elements,  and  [u ]qK+k 
is  the  measured  pressure  at  time  t  corresponding  to  the  kth  temporal  sample  recorded  by  the  qth  transducer 
element.  In  Eqn.  (1),  de(t)  is  the  time  derivative  of  the  acoustic-electric  response  (EIR)  of  the  transducer,  and  it 
is  assumed  to  be  identical  for  all  transducer  elements.  Also,  g(rs,t)  is  related  to  the  spherical  Radon  transform 
(SRT)  by  the  following  expression: 


g(rs,t) 


P  g{rs,t) 

AnCp  t 

jvdrA(v)6(c0t-\v‘  -v\), 


(2) 

(3) 


where  g(rs,£)  is  the  SRT  calculated  at  time  t  for  the  transducer  element  located  at  position  rs,  A(r)  is  the 
absorbed  optical  energy  density,  and  /3,  cq,  and  Cp  are  the  coefficient  of  volume  expansion,  the  speed  of  sound 
in  the  medium,  and  the  specific  heat  capacity  of  the  medium  at  constant  pressure,  respectively. 


2.2  Interpolation-based  Discrete-to-Discrete  (D-D)  Imaging  Model 

To  perform  iterative  image  reconstruction,  we  discretized  the  object  A( r)  in  Eqn. (3).  In  this  work,  we  have 
focused  on  the  interpolation-based  imaging  model  for  discretization.15  We  have  included  a  brief  description  of 
the  imaging  model  for  completeness  (Ref.  [15]  provides  a  more  detailed  description).  To  apply  iterative  image 
reconstruction  algorithms,  we  employed  the  following  N-dimensional  representation  of  the  object  function 


N  —  l 


-4(r)  «  y^[a]„0n(r), 

n= 0 


(4) 


where  a  is  a  coefficient  vector  whose  n-th  element  is  denoted  by  [a]n,  and  0n(r)  is  an  expansion  function.  The 
coefficient  vector  elements  are  defined  as  samples  of  the  object  function  on  the  nodes  of  a  uniform  Cartesian 
grid: 

[<*]„=[ dr 8 (r -rn)A(r),  n  =  0, 1,  ■  •  •  , N  -  1,  (5) 

Jv 

where,  rn  =  {xn,yn,zn)T  specifies  the  location  of  the  n-th  node  of  the  uniform  Cartesian  grid.  The  definition  of 
the  expansion  function  depends  on  the  choice  of  interpolation  method.17  For  a  trilinear  interpolation  method, 
the  expansion  function  can  be  expressed  as18 


lint 


(1  -  -  ^)(l  -  ^),  if  k  -  xn\,  I y  -  yn I,  \z  -  zn\  <  A. 

0,  otherwise 


(6) 


where  As  is  the  voxel  size.  On  substitution  from  Eqn.  (4)  into  Eqn.  (1),  a  D-D  mapping  from  a  to  u  can  be 
expressed  as 

u  «  Ha.  (7) 


H  can  be  decomposed  in  two  steps: 

u  =  Hqe  DeGa,  (8) 

where  G  and  De  are  discrete  approximations  of  Eqn.  (3)  and  the  operator  that  implements  a  temporal  convolution 
with  the  derivative  of  the  EIR,  respectively.  We  implemented  G  using  a  ‘ray-driven’  implementation19-21  of  the 
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SRT.  In  this  case,  for  each  data  sample,  the  surface  integral  is  evaluated  as  the  sum  of  the  contributions  from 
the  voxels  that  reside  on  the  spherical  shell  surface  specified  by  the  data  sample.  By  use  of  Eqns.  (3),  (4),  (5), 
and  (6),  we  obtain 


N~1  Ni-lNj-1 

i^^qK+k  =  5-/  Mn  ^2  ^2  0n(Uc,2,j)  =  [g]g#+fc> 

t  n= 0  2—0  j= 0 


(9) 


where  [g\qK+k  ~  0(rg,£)lt=fcAt  with  r*  specifying  the  location  of  the  q- th  point-like  transducer,  and  Ni  and  Nj 
denote  the  numbers  of  divisions  over  the  two  angular  coordinates  of  the  local  spherical  coordinate  system  with 
respect  to  the  position  of  the  transducer  element  (see  Wang’s  paper15  for  more  detail).  The  convolution  operator 
in  Eqn.  (1)  is  calculated  as 


[D  &\qK+k  4 TrCp^7  —  [U]qK+k’  (^) 

where  deeRfc  contains  the  temporal  samples  of  the  derivative  of  the  EIR,  and  Jr(-)  and  Jr_1(-)  denote  the  discrete 
Fourier  and  inverse  Fourier  transforms,  respectively. 


2.3  The  Voxel-Driven  Adjoint  Operator  in  the  Interpolation-based  Imaging  Model 

The  adjoint  operator  of  H  is  defined  as  IT  =  G^De^,  where 


[Detu] 


P 


qK+k  4^(7 


[T-HH*eYHn)}}qK+k  =  gqK+k- 


In  case  of  the  ray-driven  adjoint  operator 

[gV 


Q-1K-1  Ni-lNj-1 

^  22  22  \qK+k  22  22  =  [a  ]n- 

q— 0  k= 0  2—0  j— 0 


(ii) 


(12) 


Here,  with  respect  to  the  projection  operator,  we  introduce  an  unmatched  adjoint  operator  that  accumulates 
the  contribution  from  all  the  transducers  to  a  given  voxel.  This  scheme  is  referred  to  as  the  voxel-driven  adjoint 
operator.  For  the  voxel-driven  adjoint  operator, 

Q- 1 

[Gtoxg']n  =  A^[g']-^Kox]n.  (13) 

g=0 

Here  k  =  (|r®  —  rn|/co)/At,  and  the  value  of  [g'j^  was  approximated  by  linear  interpolation  from  its  two 
neighboring  samples  as 

[s']fc  =  k  [{k  +  1  -  k)  [g']qK+k  +  (k-k)  [g']qK+k+1] ,  (14) 

Lk 

where  =  |r*  —  rn|/co,  and  k  is  the  integer  part  of  k.  Comparing  Eqn.  (12)  with  Eqn.  (13),  we  find  that  [a'vo x]n  « 
[a']n.  Moreover,  Eqn.  (13)  can  also  be  derived  from  Eqn.  (12).  In  this  study,  the  accelerated  interpolation-based 
imaging  model  refers  to  the  D-D  imaging  model  with  forward  operators  in  Eqn.  (9)  and  (10)  and  the  adjoint 
operators  described  in  Eqn.  (11)  and  (13). 

All  operators  were  implemented  using  the  CUDA  framework  to  take  advantage  of  GPU  parallelization.  The 
calculations  were  performed  primarily  on  NVIDIA  K40  GPUs. 
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3.  COMPUTER  SIMULATIONS  AND  EXPERIMENTAL  STUDIES 


Computer  simulation  studies  were  performed  on  a  3D  phantom  containing  nine  spheres  of  different  diameters  and 
absorbed  energy  densities.  A  2D  slice  of  the  phantom  is  shown  in  Fig.  1(a).  The  pressure  data  were  generated  for 
90  tomographic  views  distributed  uniformly  between  0  and  27r  for  a  transducer  array  consisting  of  64  transducer 
elements.  The  elements  were  equiangular ly  distributed  between  0  and  i r.  The  phantom  had  an  overall  size  of 
29.4  x  29.4  x  61.4  mm3,  and  we  used  a  voxel  size  of  0.14  mm  in  each  dimension. 

Experimental  studies  were  conducted  on  pre-existing  mouse  data  from  180  tomographic  views  and  63  trans¬ 
ducer  elements.  The  transducer  elements  were  distributed  uniformly  between  14°  and  84°.  The  measured  data 
consisted  of  1536  temporal  samples  recorded  at  a  20  MHz  sampling  rate.  In  this  study,  we  set  Co  to  be  1.54 
mm//is  and  the  Griineisen  coefficient  (^Co/Cp)  to  be  2.3716. 

3.1  Image  Reconstruction  Algorithm 

3.1.1  Penalized  Least-Squares  with  Quadratic  Penalty  (PLS-Q) 

The  absorbed  optical  enery  density  was  estimated  by  minimizing  the  penalized  least-squares  objective  function 
using  a  linear  conjugate  gradient  method: 

d  =  argmin  ||u  —  Ha||2  +  /?R(a),  (15) 

OL 

where  is  the  regularization  parameter  and  R(a)  is  defined  as 

N-2 

R(o:)  =  ^  ^  [(^n  Qnx)  T  iS^n  & ny )  T  {oLn  OL nz )  ]  •  (10) 

n— 0 

In  the  above  equation,  nx ,  ny,  and  nz  are  the  neighboring  grid  points  to  the  nth  node  along  each  direction. 


4.  RESULTS  AND  DISCUSSION 

The  images  reconstructed  from  the  numerical  phantom  data  are  shown  in  Fig.  1.  It  can  be  seen  that  the 
image  from  the  new  interpolation-based  algorithm  (c)  agrees  very  well  with  that  reconstructed  using  the  original 
algorithm15  (b),  as  well  as  with  the  phantom  (a).  In  this  study,  a  penalty  was  not  employed,  i.e.  f3  =  0.  The 
images  were  obtained  after  the  34th  iteration.  Using  the  original  interpolation-based  algorithm,  it  took  918 
minutes  (using  one  GPU)  to  obtain  an  image,  while  the  acceleratd  algorithm  provided  an  image  in  41  minutes 
(using  four  GPUs).  The  results  suggest  that  the  new  accelerated  algorithm  is  20  times  more  efficient  than  the 


0.75  ^ 
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—  —  Ray-driven 
. Voxel-driven 


(d) 


ru 


in 
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Figure  1.  Slices  of  (a)  the  numerical  phantom  and  the  images  reconstructed  using  (b)  the  interpolation-based  algorithm 
with  the  ray-driven  adjoint  operator  and  (c)  the  interpolation-based  algorithm  with  the  voxel-driven  adjoint  operator;  (d) 
profiles  corresponding  to  the  vertical  arrow  in  (a)  for  all  three  images;  (e)  profiles  corresponding  to  the  horizontal  arrow 
in  (a)  for  all  three  images. 
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(a) 


(b) 


Figure  2.  Maximum  intensity  projection  renderings  of  3D  images  of  a  mouse  body  using  (a)  the  original  interpolation-based 
algorithm,  and  (b)  the  accelerated  interpolation-based  algorithm. 

original  implementation  when  we  take  advantage  of  multiple  GPUs,  and  more  than  five  times  more  efficient 
when  using  a  single  GPU  implementation. 

Images  reconstructed  for  mouse  data  are  shown  in  Fig.  2.  The  regularization  parameter  was  0.1  and  the  images 
were  obtained  after  30  iterations.  By  comparing  the  images,  we  see  that  the  accelerated  algorithm  retains  the  fine 
vasculature  structures  of  the  image  reconstruced  using  matched  ray-driven  backprojection  algorithm.  It  took  26 
hours  and  30  minutes  for  the  image  reconstruction  using  the  original  algorithm,  but  the  image  was  reconstructed 
in  just  one  hour  and  21  minutes  using  the  accelerated  algorithm. 

5.  SUMMARY 

In  this  study,  we  have  demonstrated  the  use  of  an  unmatched  adjoint  operator  (with  respect  to  the  projection 
operator)  for  both  a  numerical  phantom  and  experimental  studies.  The  unmatched  adjoint  operator  is  five  times 
more  computationally  efficient  than  the  matched  adjoint  operator.  Moreover,  the  unmatched  operator  can  also 
be  more  conveniently  parallelized  in  multi-GPU  implementations.  The  study  beneficially  extends  3D  iterative 
image  reconstruction  algorithms  to  applications  with  a  large  FOV,  such  as  breast  imaging. 
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Abstract 

In  this  work,  we  introduce  a  novel  three-dimensional  imaging  system  for  in  vivo  high-resolution 
anatomical  and  functional  whole-body  visualization  of  small  animal  models  developed  for 
preclinical  and  other  type  of  biomedical  research.  The  system  (LOUIS-3DM)  combines  a 
multiwavelength  optoacoustic  tomography  (OAT)  and  laser-induced  ultrasound  tomography 
(LUT)  to  obtain  coregistered  maps  of  tissue  optical  absorption  and  speed  of  sound,  displayed 
within  the  skin  outline  of  the  studied  animal.  The  most  promising  applications  of  the  LOUIS- 
3DM  include  3D  angiography,  cancer  research,  and  longitudinal  studies  of  biological  distributions 
of  optoacoustic  contrast  agents. 

Keywords 

photoacoustic  tomography,  speed  of  sound  tomography,  mouse  imaging 


Introduction 

Three-dimensional  full-aperture  optoacoustic,  also  known  as  photoacoustic,  tomography  (3D- 
OAT)  was  previously  shown  to  efficiently  visualize  central  and  peripheral  vasculature  and  blood- 
rich  organs  of  a  mouse.1-3  Reconstructed  3D-OAT  images  could  provide  information  valuable  for 
preclinical  studies  of  pathologies  that  affect  local  blood  content,  such  as  cancer,  trauma,  isch¬ 
emia,  stroke,  and  so  on.4’5  For  example,  in  cancer  research  development  of  angiogenesis,  a  com¬ 
plex  microvascular  network  feeding  a  growing  tumor  can  be  assessed  with  3D-OAT.4  Other 
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c  =  1.522  mm /(.is  c  =  1.526  mm/|Lis 


Figure  I.  Section  of  a  3D  optoacoustic  image  of  two  hairs.  Data  were  acquired  by  using  a  64-channel 
arc-probe,  25MHz  sampling  rate,  and  150  angular  acquisitions.  Hairs  were  aligned  in  a  cross-pattern 
stretched  along  the  axis  of  rotation,  (a)  Image  reconstructed  with  correct  speed  of  sound  of  1 .522  mm/ 
ps,  corresponding  to  the  speed  of  sound  in  the  ambient  distilled  water  at  36°C.  Image  shows  sharp 
cross-sections  of  both  hairs  with  diameters  of  260  pm.  (b)  Image  reconstructed  from  the  same  data  with 
speed  of  sound  increased  by  0.3%  to  1.526  mm/ps.  Blurring  of  the  hairs  is  noticeable,  and  the  measured 
FWHM  diameter  of  both  hairs  is  increased  by  140%  to  630  pm. 


preclinical  applications  of  full-aperture  optoacoustic  systems  include  the  monitoring  of  high- 
frequency  ultrasound  (Hi-Fu)  treatment,  longitudinal  studies  of  postmortem  changes,  and  mea¬ 
surements  of  biodistributions  of  optoacoustic  contrast  agents  such  as  carbon  nanotubes,  metal 
plasmonic  nanoparticles,  and  fluorescent  proteins.6-13  Clinical  applications  of  3D-OAT  are  also 
being  investigated  and  may  require  high-quality  visualization  of  peripheral  vasculature  down  to 
the  microcirculation  level.  One  such  example  is  the  monitoring  of  an  individual’s  peripheral 
blood  flow  during  a  physiological  stress  test  aimed  to  study  local  vasomotor  response  during 
induced  short-term  hypothermia.14  The  feasibility  of  using  3D  optoacoustic  tomography  (OAT) 
for  breast  angiography  was  also  reported.15 

Biological  objects  visualized  with  OAT  typically  consist  of  blood  vessels,  muscles,  connec¬ 
tive  tissue,  bones,  air,  and  liquid  filled  spaces.  Such  mechanically  heterogeneous  composition 
complicates  optoacoustic  imaging  by  creating  artifacts  due  to  mismatch  in  acoustic  impedances.16 
It  also  decreases  image  quality  through  blurring  introduced  by  assumption  of  homogeneous 
speed  of  sound  (SOS)  during  tomographic  reconstruction.16’17  Mitigation  of  OAT  image  artifacts 
due  to  unmodeled  SOS  variations  is  particularly  important  in  the  case  of  3D-OAT.  In  such  cases, 
the  reconstructed  objects  can  be  located  several  centimeters  from  the  acoustic  transducers,  and 
variations  in  the  SOS  introduce  noticeable  blurring  on  the  reconstructed  image.  Figure  1  shows 
an  example  of  a  section  across  3D  image  of  two  hairs  reconstructed  with  two  SOSs  different  by 
just  0.3%.  In  this  situation,  small  blood  vessels  can  be  visually  blurred  and  their  diameters 
increased  by  140%  to  630  pm  from  the  original  260  pm. 

Ultrasound  computed  tomography  (UCT)  seeks  to  reconstruct  distributions  that  describe 
acoustic  properties  of  an  object.18’19  In  biomedical  applications,  UCT  has  been  employed  to  esti¬ 
mate  the  SOS  and  acoustic  attenuation  (AA)  distributions,  which  can  characterize  soft  tissues. 
Most  previous  investigations  of  UCT  have  utilized  electrically  generated  ultrasound.20-22  UCT 
continues  to  mature  and  is  currently  being  evaluated  for  breast  cancer  imaging.  For  example,  Li 
et  al.  (2009)  used  the  CURE  system  with  256  ultrasound  transducers  evenly  distributed  over  a  20 
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Figure  2.  (a)  A  schematic  showing  a  device  for  generation  of  planar  monopolar  laser  ultrasound  waves 
(side  view).  (I)  optically  absorbing  layer  (black  paint),  (2)  acoustically  matching  backing,  (3)  optical  fiber 
bundle,  (4)  lens,  (5)  expanding  laser  beam.  Computer  simulation  of  ultrawide  band  laser  ultrasound 
planar  waveforms  generated  by  a  1 6  ns  laser  pulse  incident  upon  a  highly  absorbing  flat  substrate 
with  acoustically  matching  conditions:  (b)  temporal  profiles  and  (c)  amplitude  spectra  of  laser-induced 
ultrasound  signals  (OA  signal)  measured  at  two  different  distances  (z)  from  the  laser  ultrasound  emitter. 
D  is  a  diameter  of  the  Gaussian  laser  beam  at  the  surface  of  the  emitter. 


cm  diameter  ring  array  to  reconstruct  parallel  coronal  SOS  slices  in  61  breast  cancer  patients.23 
Their  advanced  tomography  algorithm  based  on  the  iterative  bent-ray  model  regularized  with 
total-variation  (TV)  allowed  imaging  and  quantification  of  SOS  in  malignant  lesions  (1548  ±  17 
m/s),  as  well  as  differentiation  between  fatty  tissue  (1422  ±  9  m/s)  and  breast  parenchyma 

(1487  ±  21  m/s).  They  also  demonstrated  that  SOS  images  fused  with  back-scattered  ultrasound 
could  be  used  for  monitoring  clinical  response  of  breast  cancer  to  neo-adjuvant  chemotherapy. 
Nebeker  and  Nelson  (2012)  proposed  another  technological  approach  to  coronal  breast  UCT  by 
positioning  a  breast  between  a  regular  convex  ultrasound  transducer  and  a  concave  ultrasound 
reflector,  and  engaging  a  rotational  mechanism  to  acquire  the  required  set  of  views.24 

Laser  ultrasound  (LU)  pulses  are  generated  via  the  same  optoacoustic  effect  that  allows  opto- 
acoustic  imaging.4’25  The  main  difference  is  that,  in  case  of  laser-induced  ultrasound,  absorption  of 
optical  energy  occurs  in  an  artificial  emitter  element,  which  is  located  outside  of  the  interrogated 
tissue  or  object  and  has  extremely  high  absorption  coefficient  (103  cm-1  and  more)  and  increased 
efficiency  of  thermoacoustic  conversion  (Gruneisen  parameter).  One  of  the  easiest  ways  to  create 
an  efficient  laser  ultrasound  source  is  to  use  a  Q-switched  laser  (same  laser  could  be  used  in  opto¬ 
acoustic  imaging)  illuminating  black-colored  or  painted  polymers,  such  as  polyethylene,  acrylic, 
or  Teflon,  which  have  Gruneisen  parameter  two  to  four  times  that  of  water.  High-quality  ultrawide 
band  nonreverberating  acoustic  pulses  are  achieved  by  acoustic  matching  and/or  removal  of  the 
emitter’s  interfaces.26-28  Figure  2  shows  results  of  a  computer  simulation  for  ultrawide  band  laser 
ultrasound  waveforms  generated  by  a  16  ns  laser  pulse  incident  upon  a  highly  absorbing  flat  sub¬ 
strate  with  acoustically  matching  conditions.  As  it  is  seen  in  Figure  2(b)  and  (c),  when  properly 
designed,  laser  ultrasound  sources  can  generate  acoustic  pulses,  which  are  by  far  superior  to  those 
of  traditional  electrically  generated  ultrasound  in  terms  of  lack  of  reverberations  and  spectral 
bandwidth  providing  the  best  source  for  SOS  and  spectral  AA  tomographies.27’29  In  addition,  using 
laser-induced  ultrasound  is  the  most  efficient  and  inexpensive  way  to  combine  UCT  and  OAT,  as 
both  technologies  employ  the  same  optimized  lasers,  transducers,  and  electronics. 

Two-dimensional  laser  ultrasound  tomography  (LUT)  was  originally  proposed  by  Manohar 
et  al.  (2007). 30  By  introducing  a  graphite  rod  into  the  optical  beam  of  their  transmission-mode  2D 
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photoacoustic  tomography  system  (between  the  1064  nm  laser  illuminator  and  the  interrogated 
object),  the  authors  were  able  to  separate  in  time  photoacoustic  transients  and  laser  ultrasound 
pulses  generated  by  the  illuminated  carbon  fiber.  SOS  images  of  an  agar/oil  phantom  were  recon¬ 
structed  using  straight-ray  fan-beam  CT  and  showed  good  quality  of  two  8  mm  and  5  mm  square 
inclusions.  Subsequent  publications  from  the  same  research  group  further  advanced  the  technique 
by,  first,  implementing  the  concomitant  photoacoustic  tomography  and  LUT  for  both  SOS  and 
AA,31  and,  then,  by  significantly  accelerating  data  acquisition  via  multiple  (nine)  simultaneously 
emitting  LU  elements.32  The  latter  case,  however,  resulted  in  image  degradation  as  compared  with 
a  unit  armed  with  just  a  single  LU  emitter.  The  most  recent  work  in  the  field  of  2D-LUT  demon¬ 
strated  coregistered  photoacoustic  and  SOS  maps  from  the  excised  kidney  of  a  mouse  embedded 
in  an  agar  matrix.17  Although,  high  fidelity  photoacoustic  images  showed  internal  structure  of  the 
kidney,  SOS  map  was  only  able  to  segment  the  entire  organ  from  the  background. 

In  case  when  LUT  is  coregistered  with  optoacoustic  data,  the  spatial  distribution  of  acoustic 
properties  can  be  employed  in  iterative  reconstruction  algorithms  to  improve  the  fidelity  of  opto¬ 
acoustic  images.33-35  Recently,  Xia  et  al.  (2013)  used  experimental  data  obtained  from  a  mouse 
kidney  phantom  and  from  a  leaf  phantom  to  compare  photoacoustic  reconstructions  done  with  a 
single  SOS  value  optimized  for  the  entire  interrogated  volume  with  those  utilizing  the  estimated 
map  of  SOS  distribution.17  The  authors  concluded  that  the  former,  simplified  approach  resulted 
in  images  of  reasonable  quality,  providing  spatial  SOS  variations  were  not  too  large  (as  in  the 
mouse  kidney  phantom).  When  spatial  SOS  variations  were  significant  (leaf  phantom),  photo¬ 
acoustic  images  reconstructed  with  optimized  single  SOS  demonstrated  inferior  quality,  mostly 
at  the  periphery  of  the  rotational  reconstruction  geometry. 

In  this  work,  we  present  novel  instrumentation  and  a  tomographic  framework  for  performing 
3D-OAT  with  coregistered  3D  laser-induced  ultrasound  tomography  (3D-LUT).  The  ultimate 
goal  of  this  study  is  to  develop  a  3D  laser  OAT-UCT  imaging  system  for  preclinical  research 
(LOUIS-3DM).  We  report  the  first  fully  3D  laser  UCT  unit  and  demonstrate  the  performance  of 
the  system  prototype  through  multiwavelength  optoacoustic  imaging  coregistered  with  SOS 
imaging  of  a  phantom. 

Materials  and  Method 

Optoacoustic  Imaging  Unit 

The  optoacoustic  imaging  component  of  the  system  was  based  on  a  design  proven  to  produce 
excellent  volumetric  images  of  mice  in  vivo.1’8  Figure  3a  depicts  the  imaging  module  that  con¬ 
tains  a  metal  housing  (1),  array  of  ultrasonic  transducers  (2),  alignment  unit  (3),  optical  fiber 
bundles  for  optoacoustic  imaging  of  the  skin  (4),  optical  fiber  bundles  for  optoacoustic  imaging 
of  deep  tissue  (5),  array  of  laser  ultrasound  emitters  (6),  and  computer  controlled  rotational  stage 
(7).  In  these  studies,  we  used  a  water-tank-based  prototype  of  LOUIS-3DM  (Figure  3b).  Using  a 
rectangular  water  tank  allowed  us  to  optimize  parts  and  units  during  development  process  with¬ 
out  expensive  remanufacturing  of  the  imaging  module.  The  prototype  contained  an  arc-shaped 
array  (2)  (Imasonic  SAS,  Voray  sur  l’Ognon,  France)  of  64  (2  mm  x  2  mm)  piezo-composite 
ultrasonic  transducers  (1.5-4. 5  MHz  bandwidth  at  -6dB),  which  was  used  in  both  optoacoustic 
and  ultrasound  tomography  modes.  The  aperture  of  the  array  was  spanning  a  150°  arc  with  a 
radius  of  65  mm,  oriented  vertically  inside  the  tank.  During  the  OAT  scanning,  a  mouse  or  a 
phantom  was  rotated  about  the  axis  passing  through  the  center  of  the  array  parallel  to  the  plane 
of  the  arc.  Since  all  the  transducers  in  such  a  scanning  geometry  look  straight  into  the  center  of 
the  reconstructed  volume,  the  negative  contribution  of  the  transducer’s  directivity  into  the  quality 
of  reconstructed  images  is  reduced.  Alignment  of  the  axis  of  rotation  was  performed  as  needed 
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Figure  3.  (a)  Designer’s  vision  of  the  3D-OAT/LUT  imaging  module,  (b)  Water-tank- based  prototype 
of  the  LOUIS-3DM.  Metal  housing  (I),  array  of  ultrasonic  transducers  (receivers)  (2),  manual 
alignment  unit  (3),  optical  fiber  bundles  for  optoacoustic  imaging  of  the  skin  (4),  optical  fiber  bundles 
for  optoacoustic  deep  tissue  imaging  (5),  array  of  laser  ultrasound  emitters  (6),  computer  controlled 
rotational  stage  (7),  motorized  linear  stage  (8),  gas  anesthesia  delivery  module  (9).  OAT  =  optoacoustic 
tomography;  LUT  =  laser-induced  ultrasound  tomography. 


using  a  manual  2D  translation  stage  (3)  (Thor  Labs,  Newton,  New  Jersey,  USA)  oriented  perpen¬ 
dicularly  to  the  arc  of  the  array.  The  alignment  was  confirmed  using  optoacoustic  imaging  of  a 
horse  hair. 

The  OAT  scan  was  performed  by  a  DC  motor  with  an  optical  encoder  (7)  (Faulhaber  Gmbh, 
Schoenaich,  Germany)  for  10  full  360°  rotations  or  cycles  with  2.4°  steps  and  150  views  per 
cycle.  The  rotational  step  size  was  set  to  match  the  angular  pitch  between  the  elements  of  the  arc 
array.  A  motorized  linear  stage  (8)  provided  up  to  150  mm  movement  along  the  axis  of  rotation 
and  served  simply  to  position  the  interrogated  sample  at  an  appropriate  depth  inside  the  tank.  The 
mouse/phantom  holder  with  the  gas  anesthesia  delivery  module  (9)  (Summit  Anesthesia 
Solutions,  Bend,  Oregon)  was  described  in  detail  in  our  previous  work.1  Optoacoustic  illumina¬ 
tion  was  provided  by  two  bifurcated,  randomized  fiber  bundles  with  four  rectangular  output 
profiles  of  1  mm  x  80  mm  each.  A  Q-switched  laser  system  (Spectra Wave,  Tomowave 
Laboratories,  Houston,  Texas,  USA)  was  used  to  establish  appropriate  optoacoustic  illumination 
pattern  selected  from  two  fixed  (532  nm  and  1064  nm,  9  ns  pulses,  up  to  160  mJ/pulse)  and  a 
tunable  (750-830  nm,  12  ns  pulses,  up  to  80  mJ/pulse)  output  wavelength,  and  operated  at  10  Hz 
pulse  repetition  frequency.  A  532  nm  wavelength,  being  absorbed  primarily  within  a  narrow 
superficial  layer  of  skin,  was  used  in  a  backward  optoacoustic  mode  (4)  to  outline  the  visualized 
biological  object  and  to  identify  two  separate  acoustic  domains  (inside  and  outside  of  the  biologi¬ 
cal  object)  for  future  improved  implementation  of  OAT/LUT.1617  The  wavelengths  of  760  nm, 
800  nm,  and  1064  nm  are  typically  used  in  orthogonal  optoacoustic  mode  (5),  as  they  are  neces¬ 
sary  for  reconstruction  of  functional  volumetric  maps  related  to  distribution  of  total  hemoglobin, 
water,  and  local  levels  of  blood  oxygenation.36’37  Electronic  system  components  like  40  to  90  dB 
broadband  analog  amplifier  with  adjustable  gain  control  and  the  digital  acquisition  (DAQ)  hard¬ 
ware  operated  at  64  parallel  channels,  25  MHz,  and  1536  samples  per  channel  were  previously 
reported  for  3D-OAT  mouse  imaging  system.1’8 
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Signal  conditioning  for  deep  tissue  imaging  included  (a)  synchronization  with  the  moment  of 
laser  emission  using  laser  pulse  signals  recorded  by  a  photodiode,  (b)  averaging  data  acquired  at 
identical  rotational  positions  during  the  scan,  (c)  Wiener  deconvolution  of  the  system’s  acousto¬ 
electric  impulse  response,  and  (d)  bandpass  filtering  with  0.25  to  5  MHz  bandwidth.  Data  col¬ 
lected  after  each  complete  deep  tissue  OAT  scan  were  used  in  near  full-aperture  3D  filtered 
back-projection38  for  reconstruction  of  rectangular  volume  with  0.1  mm3  voxels  and  dimensions 
up  to  30  mm  x  30  mm  x  50  mm  centered  at  the  focus  of  the  probe.  It  was  developed  to  run  on 
graphics  processing  units  (GPUs)  which  support  compute  unified  device  architecture  (CUD A), 
resulting  in  fast  reconstruction  of  3D-OAT  images  (about  1.5  min  for  50  million  voxels).39 
Visualization  was  performed  using  VolView  2.0  (Kitware,  Clifton  Park,  New  York,  USA)  and 
included  (a)  median  filtering  with  a  3  x  3  x  3  voxel  kernel,  and  (b)  two  linear  ramps  for  the  Scalar 
Opacity  Mapping  starting  near  the  mode  of  the  histogram.  Gradient-based  opacity  was  manipu¬ 
lated  manually  to  clean  out  fuzziness/noise  and  emphasize  large  intensity  changes  (boundaries). 

Signal  conditioning  for  skin  imaging  included  (a)  enveloping,  (b)  low-pass  filtering,  and  (c) 
limiting  data  to  the  first  and  last  sample  of  the  mouse  (skin)  signal  using  noise-based  threshold¬ 
ing.  The  OAT  was  then  performed  in  the  same  manner  as  in  the  case  of  deep  tissue  imaging. 
Visualization  of  the  skin  also  was  a  double  ramp  process.  Skin  volume  was  presented  in  gray 
palette,  with  a  dark  gray  scale  of  0.3  starting  near  the  mode  of  the  histogram.  Gradient  opacity 
was  manually  set  in  attempt  to  enhance  the  skin  layer. 

For  combined  representation  of  the  skin  over  internal  structures  of  a  mouse,  the  volumes  of 
skin  (532  nm)  and  blood  vessels  (1064  nm)  were  merged  together  in  Volview.  They  were  visually 
processed  in  the  same  manner  as  they  were  individually  with  the  only  exception  that  the  532  nm 
volume  was  made  70%  more  transparent  to  see  the  1064  nm  volume  within.  The  skin  was  repre¬ 
sented  in  gray  scale  palette  and  essentially  enveloped  the  1064  nm  yellow/red  image  providing 
natural  anatomical  fiducials  for  vascular  network,  kidneys,  and  spleen. 

Laser-Induced  Ultrasound  Imaging  Unit 

Several  designs  of  the  array  of  laser  ultrasound  emitters  were  considered,  and  parameters  were 
optimized  through  both  computer  simulations  and  prototyping,  while  the  design  of  the  arc  array 
of  acoustic  receivers  (the  same  that  is  used  in  the  optoacoustic  mode)  was  fixed.  The  recon¬ 
structed  volume  was  constrained  to  be  equal  to  30  mm  x  30  mm  x  50  mm  as  in  typical  whole- 
body  3D-OAT  of  a  mouse. 

As  our  computer  simulations  suggested  that  SOS  images  were  comparable  for  the  spiral,  oblique 
circular,  and  linear  diagonal  arrangements  of  LU  emitters,40  we  decided  to  further  pursue  linear 
geometry  that  could  be  accurately  reproduced  and  aligned.  Section  SI  in  the  Supplementary 
Materials  contains  details  of  the  design  and  evaluation  of  another  promising  prototype  based  on  the 
spiral  arrangement  of  the  LU  emitters.  We  fabricated  a  rectangular  acrylic  support  of  optical  fibers 
with  fiber  slots  arranged  in  a  planar  diagonal  pattern  facing  the  array  of  receivers  (Figure  4a).  The 
output  termini  of  thirty-three  600  pm  multimode  step-index  fibers  were  permanently  glued  to  the 
acrylic  support.  The  inputs  of  the  fibers  were  arranged  in  a  linear  pattern  and  fixed  within  the 
aluminum  plate  (Figure  4b).  Both  ends  of  the  fiber  were  polished  to  0.3  micron  with  lapping  film. 
During  LUT  scans,  fiber-optic  inputs  of  the  emitters  were  coupled  to  the  fiber-optic  output  of  the 
laser  using  a  ball-lens  projection  system  (Figure  4c).  A  cheap  and  disposable  laser  printer  trans¬ 
parency  film  (thickness  about  125  pm)  painted  with  black  plastic-specific  spray  acrylic  was  used 
to  produce  laser  ultrasound  pulses.  The  film  was  positioned  53  mm  from  the  axis  of  rotation  with 
painted  side  facing  the  array  of  receivers.  To  minimize  acoustic  reflections  generated  on  a  glass 
interface,  the  painted  film  was  offset  by  8  mm  with  respect  to  the  fiber-optic  illuminators  using 
two  rectangular  acrylic  spacers.  Due  to  the  short  offset  distance  and  small  difference  between 
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Figure  4.  Laser-induced  ultrasound  unit,  (a)  Array  of  LU  emitters,  (b)  Fiber-to-fiber  optic  coupling  and 
scanning,  (c)  Ball-lens  (LMS-LSFN)  fiber-optic  coupler.  (I)  Acrylic  support  of  optical  fibers,  (2)  thirty- 
three  600  pm  optical  fibers,  (3)  transparency  film  painted  with  black  acrylic  (LU  source),  (4)  aluminum 
support  of  fiber-optic  inputs  of  the  emitters,  (5)  fiber-optic  output  of  the  laser,  (6)  ball-lens  projection 
system,  (7)  linear  motorized  stage. 


refractive  indices  of  the  fiber’s  core  and  water,  light  beam  expansion  was  negligible.  Each  indi¬ 
vidual  LU  emitter  mimicked  a  point  source  by  exploiting  the  expansion  of  a  very  small  diameter 
(600  pm)  flat  wave.  We  used  fundamental  Nd:YAG  wavelength  of  the  Spectra  Wave  laser  system 
attenuated  to  energies  about  100  pJ  per  pulse.  Laser  ultrasound  was  generated  by  conversion  of 
9  ns  laser  pulses  into  clean  nonreverberating  ultrawide  band  acoustic  pulses. 

The  temperature  of  water  in  the  tank  (background)  was  maintained  by  heating  elements  placed 
at  the  bottom  of  the  tank  and  controlled  by  the  proportional-integral-derivative  controller  (PID 
Temperature  Controller,  Watlow  Inc.,  Columbia,  Missouri,  USA).  The  precision  of  ±0.1  °C  was 
achieved  during  both  calibration  and  experiment  facilitating  accurate  estimates  of  the  SOS  for 
the  background  (Figure  5).  Spatial  gradients  of  the  water  temperature  and  water  contamination 
were  minimized  by  a  custom-made  water  circulation  and  filtration  system  based  on  FilStar  XP1 
(Rena,  Chalfont,  Pennsylvania,  USA). 

At  every  rotational  position,  optical  fiber  inputs  of  laser  ultrasound  emitters  were  sequentially 
scanned  with  a  fiber-optic  output  of  the  laser  mounted  on  a  linear  motorized  stage  (Figure  4b). 
LUT  scan  was  performed  immediately  prior  to  the  optoacoustic  scans  using  deionized  water 
maintained  at  the  desired  constant  temperature  as  a  coupling  background  acoustic  medium  with 
known  speed  of  sound.  During  each  LU  emission,  all  64  receivers  acquired  1536  samples  in 
parallel  at  40  MHz  sampling  rate.  It  took  approximately  1 .5  hr  to  complete  a  scan  using  33  emit¬ 
ters,  72  rotational  views,  and  64  averaged  acquisitions/(emitter  x  view). 

Characterization  of  LU  Emitters 

A  single  element  prototype  of  the  LU  source  was  designed  and  assembled.  To  evaluate  the  ampli¬ 
tude,  spectral,  and  directivity  characteristics  of  laser  ultrasound  pulses,  the  source  was  held  at  the 
bottom  of  our  custom-built  protractor27  such  that  the  emitting  surface  was  located  exactly  at  the 
axis  of  rotation  (Figure  6a).  The  protractor  (Figure  6b)  allowed  accurate  and  reproducible  rota¬ 
tion  of  the  emitter  with  respect  to  the  axis  of  detection  of  the  calibrated  hydrophone  (Onda 
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Figure  5.  Temperature  readings  inside  the  water  tank  during  LUT  of  the  phantom.  LUT  =  laser-induced 
ultrasound  tomography. 


Figure  6.  Setup  for  characterization  of  LU  emitter,  (a)  A  prototype  of  a  single  emitter  (I)  is  mounted 
on  the  protractor  (2).  The  surface  of  the  emitter  is  located  exactly  on  the  axis  of  rotation.  The 
calibrated  hydrophone  (3)  is  used  to  measure  a  generated  acoustic  waveform,  (b)  The  protractor 
assembly.  See  text  for  further  details,  as  well  as  Ref.  27. 


HGL-1000).  The  hydrophone  had  a  broadband  (0-30  MHz)  sensitivity  of  510  nV/Pa.  The  output 
of  the  hydrophone  was  fed  to  the  oscilloscope  (Tektronix  TDS-3014)  through  the  AH-2010 
broadband  amplifier  with  a  gain  A  =  10.  The  impulses  were  acquired  with  gigahertz  sampling 
rate  (10,000  points  per  trace,  thus  10  ps  were  acquired)  and  saved  to  a  computer  with  LabView 
through  GPIB  interface.  During  acquisition,  the  emitter  and  the  hydrophone  were  separated  by 
about  50  mm  to  measure  acoustic  wavefront  in  the  far  field.  The  angle  of  detection  varied  from 
90°  (normal  to  the  surface  of  the  source)  to  40°.  To  quantify  the  efficiency  of  the  LU  source,  a 
single  point  measurement  was  also  performed  at  the  surface  of  the  emitter. 
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3D  Speed  of  Sound  Tomography 

The  goal  of  SOS  tomography  is  to  reconstruct  relative  slowness  of  the  object  with  respect  to  the 
homogeneous  medium  using  the  measured  relative  time  of  flight  (TOF)  for  all  source-receiver  pairs. 


Time  of  flight  (TOF)  estimation.  The  TOF  was  measured  using  a  steep  onset  of  a  transmitted  laser 
ultrasound  pulse  detected  within  a  time  window  defined  by  the  known  emitter-receiver  distance 
and  valid  SOS  range  (1.45-1.55  mm/ps,  considering  that  in  our  geometry  each  LU  ray  traversed 
mostly  through  ambient  water).  LU  signals  processed  with  a  bandpass  1  to  5  MHz  second  order 
Butterworth  filter  were  subject  to  thresholding  based  on  the  estimated  root  mean  square  (RMS) 
noise  (a).Matlab  filter  function  with  a  phase  recovery  was  used  to  preserve  the  steep  onset  of  LU 
signal.  The  first  sample  within  the  time  window,  which  exceeded  the  empirically  selected  thresh¬ 
old  Th  =  10a  was  selected  as  the  onset  of  the  transmitted  LU  pulse.  If  no  sample  was  found  above 
the  threshold,  the  error  was  recorded  for  that  particular  view  and  receiver-emitter  pair,  and  it  was 
not  used  in  the  subsequent  reconstruction. 

Following  mechanical  alignment  of  the  LUT  unit,  the  relative  coordinates  of  emitters  and 
receivers  were  calibrated  by  model-fitting  the  TOF  values  measured  for  background  acoustic 
medium  with  well-known  speed  of  sound.  Using  preliminary  calibration  of  TOF  for  thermally 
stable  background  medium  was  supposed  to  make  reconstruction  less  sensitive  to  system  align¬ 
ment  and  knowledge  of  absolute  coordinates  for  emitters  and  receivers. 


3D  speed  of  sound  image  reconstruction.  We  performed  image  reconstruction  based  on  the  geo¬ 
metrical  ray  theory  for  sound  waves  to  establish  a  nonlinear  model  that  relates  the  measured  TOF 
values  to  3D  SOS  distribution.  This  results  in  the  Eikonal  equation: 


VT°(F) 


(1) 


In  Equation  (1),  T°(r)  is  the  TOF  at  position  r  with  respect  to  the  emitter  location,  and 
s°(r)  is  the  slowness  distribution  and  is  the  inverse  of  the  SOS  distribution.  Our  reconstruction 
model  assumes  straight  acoustic  rays  propagating  in  3D  space.35  Based  on  the  straight-ray  model, 
for  a  given  emitter-transducer  pair  ( q ),  the  TOF,  Tq  ,  is  calculated  as  an  integral  over  the  slow¬ 
ness  distribution  along  a  straight  line: 

[T°\=\S°^)d7-  (2) 


For  a  discretized  slowness  object  model  of  dimension,  TV,  the  mapping  equation  from  object 
space  to  the  data  space  is  expressed  as 


T 


o 


q 


N- 1 


n- 0 


(3) 


where  lnq  is  the  length  segment  of  the  linear  ray  for  the  qth  emitter-transducer  pair  in  the  nth 
voxel.  The  relative  TOF  was  calculated  as  the  difference  between  the  TOF  for  a  homogeneous 
medium  of  slowness  sh  and  for  the  given  slowness  distribution.  Mathematically, 

AT  =  Th -T°  =L(sh -s°)  =  LAs,  (4) 

where  fh  is  the  TOF  vector  for  the  homogeneous  medium,  and  L  is  the  system  matrix  with 
elements  lnq •  An  estimate  of  the  relative  slowness  (a?  tf  ^  -s°)  was  reconstructed  by  solving 
the  optimization  problem, 
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Table  I.  Optical  and  Acoustic  Properties  of  the  Optoacoustic/Ultrasonic  Mouse  Phantom. 


Object 

pa(760  nm)  [cm  '] 

Ma(  1064  nm)  [cm'1] 

c  [mm/ps] 

Target  1 

2 

1 

1.42 

Target  2 

2 

0 

1.42 

Phantom  background 

0 

0 

1.57 

Acoustic  coupling  medium 

0 

0 

1.50 

As  =  min  Ills  -  Art  +  2XTV  ( As  ) , 

AseC 11  11  V  7 


(5) 


where  X  is  the  regularization  parameter.  Slowness  was  constrained  to  have  values  in  a  range 
defined  based  on  the  a  priori  knowledge  of  the  object.  We  included  the  total  variation  (TV)  regu¬ 
larization  to  accommodate  severe  incompleteness  of  the  measured  data  and  minimized  our  objec¬ 
tive  function  using  the  monotone  form  of  fast  iterative  shrinkage/threshold  algorithm  (FISTA). 
We  followed  the  implementation  of  FISTA  as  described  by  Beck  and  Teboulle.41 

Mouse 

To  demonstrate  in  vivo  multiwavelength  3D-OAT,  we  imaged  an  Athymic  Nude-Foxnlnu  mouse 
(Harlan,  Indianapolis,  Indiana)  that  was  8  weeks  old  and  weighed  approximately  30  g.  Animal 
handling,  isoflurane  anesthesia,  and  euthanasia  were  described  in  detail  in  our  previous 
publications.1’8  All  the  mouse-related  procedures  were  in  compliance  with  our  Institutional 
Animal  Care  and  Use  Committee  (IACUC)  protocol.  The  mouse  was  sequentially  scanned  using 
760  nm  for  imaging  deoxygenated  blood,  1064  nm  for  imaging  oxygenated  blood,  and,  finally, 
532  nm  for  skin  outline  imaging.  During  the  entire  scanning  procedure,  the  mouse  remained  in  the 
same  position  inside  the  holder  bracket.  The  laser  fluence  measured  at  the  animal’s  skin  was  0.4 
mJ/(pulse-cm2)  for  760  nm,  0.8  mJ/(pulse-cm2)  for  1064  nm,  and  0.5  mJ/(pulse-cm2)  for  532  nm. 

Phantom 

We  prepared  a  dual-modality  cylindrical  gelatin  phantom  with  dimensions  similar  to  those  of  a 
mouse  (025  mm)  containing  two  7  mm  spherical  polyvinyl  chloride  plastisol  (PVCP)  targets. 
Table  1  shows  measured  optical  and  acoustic  properties  of  the  phantom  and  the  aqueous  back¬ 
ground.  3D-LUT  scan  of  the  phantom  was  followed  by  two  optoacoustic  scans  at  760  nm  and 
1064  nm  laser  wavelengths  using  orthogonal  (deep  tissue)  illumination  mode. 


Results 

Performance  of  the  LUT  Unit 

The  primary  novelty  of  our  system  is  a  3D-LUT  unit,  which  was  designed  for  reconstruction  of 
volumetric  maps  of  SOS  and  spectral  AA  providing  a  new  level  of  diagnostic  information  coreg¬ 
istered  with  optoacoustic  data.  In  addition,  the  information  contained  in  ultrasonic  images  can  be 
employed  in  iterative  algorithms  of  OAT  to  improve  image  fidelity.  Laser  ultrasound  is  generated 
through  conversion  of  low-energy  (about  100  pJ)  9  ns  laser  pulses  into  clean  wideband  acoustic 
pulses  (Figure  7).  The  materials  currently  used  in  LOUIS-3DM  for  generation  of  laser  ultrasound 
provide  about  5  kPa  of  acoustic  pressure  per  1  mJ/cm2  of  incident  1064  nm  laser  radiation,  as 
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Figure  7.  A  pulse  of  laser  ultrasound  generated  by  the  LUT  prototype  emitter  (a)  and  its  spectrum 
(b)  showing  central  frequency  at  8.5  MHz  and  150%  bandwidth.  LUT  =  laser-induced  ultrasound 
tomography. 


Figure  8.  (a)  Emission  directivity  of  a  single  LU  prototype.  Graph  made  in  the  polar  coordinate  system 
shows  the  amplitude  of  the  measured  LU  signal  as  a  function  of  emission  angle  (deg).  The  data  are 
normalized  to  the  amplitude  of  the  LU  signal  measured  for  normally  emitted  LU  wave,  (b)  Variation  of 
detected  LU  signal  across  the  detector  array.  LUT  unit  provides  a  wavefront  SNR  of  5  or  greater  across 
the  entire  array.  SNR  =  signal-to-noise  ratios;  LUT  =  laser-induced  ultrasound  tomography. 


measured  with  our  calibrated  hydrophone.  Pressure  pulses  have  central  frequency  of  8.5  MHz 
and  150%  bandwidth  (Figure  7b). 

We  also  characterized  the  directivity  of  the  expanding  acoustic  waves  generated  by  our  LU 
sources.  A  true  spherical  wave  would  provide  the  best  acoustic  stimulus  for  our  LUT  system,  as 
it  would  be  detected  uniformly  along  the  entire  array  of  receivers.  Figure  8a  shows  directivity 
measured  in  the  far  field  of  the  LU  emitter,  confirming  that  unattenuated  laser  ultrasound  pulses 
(based  on  -6  dB  level)  are  generated  within  the  emission  cone  of  90°.  To  complete  the  analysis, 
emitter  performance  was  evaluated  within  normal  operating  conditions  of  the  LUT  system. 
Signals  were  acquired  for  all  detectors  with  the  central  emitter  of  the  LU  array  being  illuminated 
as  a  means  of  measuring  the  azimuthal  spread  of  the  acoustic  field.  The  impulses  were  analyzed, 
and  the  signal-to-noise  ratios  (SNR)  were  tabulated  for  all  channels.  The  data  are  plotted  as  the 
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bar  graph  shown  in  Figure  8b.  These  data  show  the  LU  emitter  excites  the  detector  array  with 
acoustic  signals  that  have  SNR  between  5  and  45,  which  is  sufficient  for  reliable  measurements 
ofTOF. 

We  performed  a  numerical  simulation  of  our  experimental  LUT  setup  and  scanning  procedure 
using  a  digital  phantom  with  3D  SOS  distribution  replicating  that  of  the  phantom  (see  the  section 
“Phantom”).  Figure  9  shows  the  results  of  the  simulation  for  the  noiseless  TOF  data — left  panel — 
and  with  additional  Gaussian  noise  (Full  width  at  half  maximum,  FWHM  =100  ns) — right  panel. 
The  images  of  the  reconstructed  meridional  slice  through  the  phantom  (bottom  panels  in  Figure  9) 
show  that  our  LUT  system  should  be  able  to  provide  high  contrast  accurate  SOS  images  even 
under  experimentally  relevant  noise  levels.  The  reconstructed  values  of  SOS  were  1 .430  ±  0.0 1 6 
mm/ps  (Average  ±  SD )  and  1.519  ±  0.009  mm/ps  for  the  spherical  inclusions  and  the  phantom 
matrix,  respectively.  The  numbers  were  in  close  agreement  with  the  original  values  of  1 .42  and 
1.52  mm/ps  used  in  the  simulations. 

Multiwavelength  3D  Optoacoustic  Imaging 

An  example  of  in  vivo  multiwavelength  3D-OAT  of  a  nude  mouse  (dorsal  view)  is  shown  in  the 
Figure  10.  The  skin  and  superficial  blood  vessels  were  visualized  using  backward  mode  illumina¬ 
tion  with  a  532  nm  laser  (Figure  10a).  Deep  tissue  imaging  using  orthogonal  illumination  with  an 
800  nm  wavelength,  which  is  close  to  the  isosbestic  point  of  hemoglobin,42  shows  details  of 
internal  vascular  network  and  blood-rich  organs,  like  kidneys,  spleen,  and  intestines  (Figure 
10b).  Using  1064  nm  laser  with  orthogonal  illumination  shifts  the  contrast  of  the  visualized  struc¬ 
tures  toward  arterial  circulation  and  tissues  with  high  water  content  (Figure  10c).  Using  our 
combined  imaging  protocol  (see  the  section  “Optoacoustic  Imaging  Unit”),  we  created  a  volume 
showing  the  internal  organs  and  circulatory  system  anatomically  associated  with  the  skin  of  the 
animal  (Figure  lOd).  A  video  showing  the  rotating  3D-OAT  image  of  the  mouse  could  be  found 
at  https://youtu.be/BzonKQqOil4. 

Coregistered  3D  Optoacoustic  and  3D  Speed  of  Sound  Imaging 

We  performed  coregistered  multiwavelength  3D-OAT  and  3D-LUT  imaging  on  a  specially 
designed  tissue  mimicking  phantom  (see  the  section  “Phantom”).  Figure  1 1  shows  a  photograph 
of  the  cylindric  gelatin  phantom  inside  the  mold  during  fabrication  process  (panel  a)  and  opto¬ 
acoustic  volumetric  maps  depicting  spherical  inclusions  according  to  their  optical  absorption 
coefficients  (panels  b  and  c).  The  reconstructed  laser-induced  ultrasound  image  is  depicted  in 
Figure  lid,  where  the  phantom  outline  and  two  inclusions  are  clearly  delineated  based  on  their 
SOS  contrast.  The  reconstructed  values  of  SOS  were  1.538  ±  0.002  mm/ps  (sphere  1),  1.546  =fc 
0.007  mm/ps  (sphere  2),  and  1.558  ±  0.018  mm/ps  (phantom  matrix).  A  video  showing  the  rotat¬ 
ing  3D  SOS  image  of  the  phantom  could  be  found  at  https://youtu.be/_PtNTN6N_VY. 

Discussion  and  Summary 

In  this  work,  we  described  the  development  of  a  3D  laser  optoacoustic  ultrasonic  imaging  system 
for  preclinical  research  (LOUIS-3DM).  A  tank-based  prototype  of  the  system  was  tested  in  phan¬ 
tom  studies  to  assess  the  quality  of  coregistered  multiwavelength  optoacoustic  and  SOS  images. 
The  obtained  OAT  and  SOS  images  are  fundamentally  different  from  those  previously  demon¬ 
strated  using  different  variants  of  a  2D-LUT  slicer.17’31’32  Jose  et  al.  (2011)  and  Resink  et  al.  (2011) 
were  able  to  visualize  2.6  mm  objects  inside  their  agar  phantoms,  while  Xia  et  al.  (2013)  were  able 
to  get  outlines  of  an  excised  kidney  of  a  mouse  (-5  x  12  mm).  However,  despite  reasonable  in-slice 
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Figure  9.  SOS  reconstruction  of  a  digital  phantom  that  consists  of  a  cylindric  24  mm  in  diameter 
body  (SOSbody  =  1.52  mm/ps)  with  two  spherical  inclusions,  10  mm  in  diameter  each  (SOSjnd  =  1.42 
mm/ps).  Acoustically  coupling  medium  has  SOScoup)e  =  1.50  mm/ps.  Top  panels  show  mapping  of  TOF 
measured  at  a  particular  rotational  view  for  all  emitter-receiver  pairs.  Due  to  the  rotational  symmetry 
of  the  phantom,  the  maps  look  identical  for  all  the  views.  Colorbars,  displayed  on  the  right  of  each  map 
indicate  TOF  deviation  (ps)  from  the  calibrated  values  (those  measured  for  the  same  emitter-receiver 
pairs  while  LU  propagated  through  background  only).  Bottom  panels  show  the  meridional  slice  through 
the  phantom  reconstructed  via  our  tomographic  algorithm  (see  the  section  “3D  Speed  of  Sound 
Tomography”).  The  inclusions  are  clearly  visible  with  high  contrast  on  both  the  noiseless  image  (left 
panel)  and  on  the  image  reconstructed  from  the  data  with  additive  100  ns  FWHM  Gaussian  noise  (right 
panel).  A  total  of  72  views  were  simulated.  Colorbars,  displayed  on  the  right  of  each  reconstructed 
image  indicate  SOS  (mm/ps).  SOS  =  speed  of  sound;  TOF  =  time  of  flight. 


resolution  of  the  SOS  images,  all  the  implementations  did  not  claim  any  practical  capability  in 
stacking  parallel  slices  of  SOS  and  AA  into  a  continuous  volume,  suffering  from  difficulties  of 
in-slice  focusing  and  interslice  registration.  Therefore,  in  these  studies,  we  proceeded  with  the 
fully  3D-LUT  approach,  which  finally  provided  high-quality  reconstruction  of  continuous  3D 
volumes  of  SOS  fully  coregistered  with  the  already  established  3D-OAT  imager,  as  it  could  be 
seen  in  Figure  S2  (Supplementary  Materials)  and  Figure  1 1 . 


Downloaded  from  uix.sagepub.com  by  guest  on  September  3,  2016 


14 


Ultrasonic  Imaging 


Figure  10.  In  vivo  3D  optoacoustic  image  of  a  nude  mouse  (dorsal  view),  (a)  Image  acquired  in  the 
backward  mode  with  532  nm  laser  illumination  showing  skin  and  superficial  blood  vessels,  (b)  and  (c) 
Deep  tissue  images  acquired  in  the  orthogonal  illumination  mode  with  an  800  and  1064  nm  wavelengths. 
Deep  tissue  OAT  images  uncover  fine  spatial  details  of  central  and  peripheral  circulatory  system  (I), 
intestine  (2),  right  kidney  (3),  and  spleen  (4).  (d)  A  composite  OAT  image  showing  spatial  details  of 
the  internal  organs  and  circulatory  system  referenced  to  the  skin  of  the  animal.  OAT  =  optoacoustic 
tomography. 


Figure  II.  (a)  Mouse  phantom  in  the  mold  with  Target  I  at  the  top  and  Target  2  below,  (b) 
Optoacoustic  image  acquired  at  760  nm  showing  both  targets,  (c)  Optoacoustic  image  acquired  at  1064 
nm  showing  only  the  top  target  (Target  I).  (d)  Reconstructed  speed  of  sound  volume  showing  the 
phantom’s  exterior  boundary  and  both  targets. 


One  obvious  disadvantage  of  our  3D-LUT  technique  as  compared  to  the  published  2D  imple¬ 
mentations  is  the  speed  of  data  acquisition.  Jose  et  al.  (2011)  employed  a  single  passive  LU  ele¬ 
ment,  45  views  per  scan,  and  100  averages  per  view,  which  totaled  to  about  8  min.31  LU  scans  by 
Resink  et  al.  (2011)  were  done  using  multielement  implementation  of  LU  emitters  (9  emitters) 
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with  otherwise  identical  setup,  which  allowed  them  to  reduce  time  of  the  scans  from  13.3  min  (80 
views)  down  to  1.5  min  (9  views),  although  with  some  degradation  of  the  image  quality.32  Our 
LUT  scans  took  1.5  hr  (72  views,  33  emitters,  64  averages)  followed  by  about  15  min  per  each 
additional  OAT  scan  (1500  views,  no  averages).  One  obvious  way  to  accelerate  our  LUT  scans  is 
to  reduce  number  of  averaging  required  for  an  acceptable  SNR.  To  do  this,  a  larger  laser  energy 
could  be  delivered  through  the  fiber-optics.  Currently,  to  prevent  crosstalk  between  adjacent  LU 
emitters,  we  coupled  well  the  output  of  the  laser  fiber  into  an  input  of  the  LU  emitter  fiber  using 
a  ball-lens  system  (Figure  4b  and  4c).  However,  the  input  of  the  600  pm  laser  fiber  is  passively 
inserted  into  a  6  mm  output  laser  beam,  which  results  in  a  significant  loss  of  the  coupled  laser 
energy  (at  least  100  times)  producing  an  output  spot  of  just  100  pJ.  Therefore,  increasing  the 
energy  of  individual  LU  sources  by  at  least  10  times  is  not  a  significant  problem  and  would 
require  an  additional  beam  coupling  optical  system  at  the  input  of  the  laser  fiber.  Stronger  emis¬ 
sion  of  LU  sources  will  also  relax  the  requirements  of  the  amplification,  therefore  reducing  the 
random  noise. 

Two-dimensional  LUT  prototypes  were  also  designed  to  simultaneously  acquire  LUT  and 
OAT  data  by  temporally  separating  OA  signals  generated  inside  the  interrogated  object  from  the 
artificial  LU  signals,  which  arrive  later.30-32  In  these  studies,  we  did  not  implement  OAT  and  LUT 
scans  to  be  run  simultaneously.  The  reason  was  to  first  demonstrate  in  a  clear  experiment  (sequen¬ 
tial  LUT  and  OAT  scans)  capabilities  and  advantages  of  the  fully  3D-LUT.  In  a  3D-OAT  geom¬ 
etry,  useful  signals  come  from  the  interrogated  object  during  the  time  interval,  which  may 
significantly  exceed  the  time  interval  for  useful  OAT  signals  in  2D  sheer.  Therefore,  the  LU 
emitters  must  be  separated  from  the  array  of  receivers  by  a  longer  distance  to  not  interfere  with 
the  OA  signals.  Longer  separation  between  the  array  of  emitters  and  receivers  increases  the 
requirements  for  the  laser  pulse  energy  to  maintain  the  same  SNR.  It  also  reduces  the  size  of 
interrogated  volume  unless  the  array  of  emitters  is  proportionally  extended.  We  still  believe  that 
possibility  of  simultaneous  OAT  and  LUT  scanning,  originally  proposed  by  Manohar  et  al.  (2007) 
for  2D  geometry,30  is  a  useful  approach  to  accelerate  and  simplify  data  acquisition,  and  we  will 
investigate  its  applicability  to  a  fully  3D-LUT  in  our  future  studies. 

Resolution  of  the  reconstructed  LUT  images  increases  with  higher  density  of  the  laser  ultra¬ 
sound  emitters/receivers.30’31’43  In  2D  geometries,  LUT  prototypes  benefit  from  the  established 
industrial  standards  for  high-quality  ultrasound  arrays  (receivers)  with  a  small  inter  element 
pitch.  In  case  of  3D-LUT,  both  emitters  and  receivers  must  be  tightly  arranged.  In  these  studies, 
our  array  of  emitters  had  33  LU  sources  separated  by  3.4  mm  from  each  other.  The  emitter¬ 
forming  optical  fibers  were  individually  embedded  into  acrylic  matrix.  To  increase  emitters’ 
density  in  the  future  prototypes,  individual  fiber  slots  could  be  abandoned  and  instead  all  the 
fibers  densely  packed  within  a  single  diagonal  groove.  Minimum  LU  emitter  pitch  for  the  same 
type  of  fibers  arranged  in  a  linear  pattern  could  be  in  this  case  reduced  to  1 . 1  mm,  allowing  three 
times  the  number  of  emitters  over  the  same  aperture  of  the  array.  Alternatively,  more  complex  2D 
patterns  of  the  LU  emitters  may  be  considered. 

Poor  directivity  of  emitter-receiver  pairs  reduces  SNR  and  widens  the  measured  LU  pulse.27  It 
eventually  leads  to  the  loss  of  fidelity  in  evaluated  TOF.  The  directivity  of  our  LUT  unit  (Figure  8) 
was  sufficient  to  provide  90°  -6dB  emission  cone  and  SNR  of  5-45  for  the  central  emitter  and  the 
entire  array  of  receivers.  For  the  best  directivity  performance,  each  emitter  should  be  aimed 
toward  the  bisector  with  respect  to  the  outermost  receivers.  However,  such  an  implementation  is 
technologically  challenging  and  was  not  attempted  during  our  prototyping.  Improved  directivity 
performance  could  be  also  achieved  by  increasing  LU  source  power  toward  the  outside  of  the 
emitters’  array  by  using,  for  example,  thinner  optical  fibers  (200  pm)  in  the  middle  and  thicker 
(600  pm)  at  the  ends  of  the  array.  Even  wider  directivity  of  individual  emitters  could  be  provided 
by  spherical  LU  sources.  Conjusteau  et  al.  (2013)  used  spherical  LU  emitters  made  out  of  ball 
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lenses  and  demonstrated  significantly  wider  emission  directivity  than  that  of  the  600  jam  planar 
source.44  However,  spherical  emitters  have  large  cross-section  and,  therefore,  impose  additional 
limitations  on  the  emitter  pitch. 

In  this  work,  we  used  3D-LUT  to  visualize  phantom  inclusions  with  SOS  10%  smaller  than  in 
the  phantom  matrix.  Our  simulations  (Figure  9)  showed  that  experimentally  relevant  100  ns  ran¬ 
dom  TOF  noise  does  not  present  a  problem  for  the  quality  of  reconstruction.  In  this  case,  both 
forward  model  of  ultrasound  wave  propagation  and  the  inverse  model  utilized  in  the  reconstruc¬ 
tion  algorithm  assumed  straight  acoustic  rays.  The  estimates  of  SOS  (for  targets  and  phantom 
matrix)  agreed  well  with  the  set  values  with  only  a  slightly  higher  number  (by  0.7%)  obtained  for 
the  SOS  of  the  targets.  However,  when  the  experimental  data  were  analyzed,  the  reconstructed 
values  for  the  SOS  in  spherical  targets  exceeded  the  measured  SOS  values  by  more  than  8%, 
while  the  SOS  estimates  remained  quite  accurate  for  the  phantom’s  matrix  (lower  by  just  0.8%). 
Both  targets  were  still  reconstructed  in  a  correct  negative  SOS  contrast  with  respect  to  the  phan¬ 
tom’s  matrix  and  allowed  clear  delineation  along  their  boundaries  (see  Figure  11).  The  discrep¬ 
ancy  in  the  reconstructed  SOS  could  originate  from  the  following  two  sources.  First,  the  utilized 
straight-ray  model  is  a  rough  approximation,  which  works  well  only  in  situations  where  SOS 
deviations  are  very  small.  When  this  assumption  is  violated,  the  acoustic  rays  are  getting  refracted 
on  heterogeneous  interfaces  and  produce  systematic  error  in  estimated  SOS.  From  this  point  of 
view,  the  SOS  of  targets  (1 .42  mm/ps)  likely  differed  enough  from  the  SOS  of  phantom’s  matrix 
(1.57  mm/ps)  to  violate  this  approximation.  Second,  the  size  of  targets  was  quite  small,  and  our 
setup  could  only  provide  very  sparse  sampling  of  the  TOF  measurements  through  the  spherical 
targets  with  few  receiver-source  pairs  providing  TOF  from  the  rays  passing  through  the  targets. 
In  this  scenario,  the  system  matrix  becomes  underdetermined,  which  we  tried  to  partially  com¬ 
pensate  by  including  a  TV  regularization.  The  penalty  term  can  also  contribute  to  the  observed 
deviations  of  SOS  from  the  expected  values. 

The  fundamental  limitations  of  TOF  sensitivity  are  defined  by  the  sampling  rate  of  the  DAQ 
system.  We  used  the  40  MHz  sampling  rate  as  Xia  et  al.  (20 13), 17  which  is  translated  into  a  maxi¬ 
mum  TOF  resolution  of  25  ns,  while  Jose  et  al.  (201 1)31  sampled  the  signals  twice  faster.  In  real¬ 
ity,  these  numbers  are  generally  superseded  by  the  accuracy  of  the  method  used  for  TOF  detection 
(see  Section  S2  in  the  Supplementary  Materials  for  more  details  on  accuracy  and  sensitivity  of 
the  TOF  measurements).  In  these  studies,  we  used  a  simple  method  based  on  the  detection  of  the 
pulse’s  onset  as  soon  as  it  exceeds  the  threshold,  which  was  established  from  the  channel’s  RMS 
noise.  Xia  et  al.  (2013)  used  cross-correlation  technique17  and  Jose  et  al.  (2011)  had  it  refined 
with  maximum  likelihood  (ML)  estimator.31  However,  both  techniques  are  still  not  ideal  in  case 
of  dispersive  behavior  of  the  acoustic  waves.  Other,  more  complex,  TOF  pickers  were  also 
developed.45  The  next  step  to  improve  accuracy  of  the  reconstructed  3D  SOS  images  would  be  to 
apply  more  sophisticated  reconstruction  algorithms.  For  example,  a  bent-ray  model  approach 
was  implemented  in  2D  SOS  tomography23  and  could  be  potentially  adapted  for  3D  problems. 

Currently,  only  3D  SOS  tomography  was  developed  in  details  and  successfully  tested  in  phan¬ 
toms.  Improved  quality  of  generated  LU  pulses  should  allow  us  in  the  future  to  make  a  next  step 
toward  reconstruction  of  spectral  ultrasound  attenuation  maps.  The  ultimate  goal  would  be  in 
vivo  coregistered  multiwavelength  optoacoustic  and  laser-induced  ultrasound  imaging  in  a 
mouse  and  modification  of  the  system  for  clinical  3D-OAT  and  UCT  applications,  like  diagnos¬ 
tics  of  breast  cancer. 
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